C.=========================================================================
C.  Library of programs for the generation of nucleus-nucleus interactions
C.  and the study of nucleus-induced cosmic ray showers
C.
C.  September 2001  changes in FPNI, and SIGMA_INI,
C.                  new SIGMA_PP, SIGMA_PPI, SIGMA_KP  (R. Engel)
C.
C.  may  1996       small bug  corrected by Dieter Heck in NUC_CONF 
C.
C.  march 1996      small modification to the superposition code
C.
C.  February 1996   change to FPNI to give an interaction length
C.                   also  at very low energy  
C.
C.  Version 1.01  september 1995 
C.       (small corrections P.L.)
C.       the random number generator is called as S_RNDM(0)
C.  ------------------------------------------------------
C.  Version 1.00  April 1992
C.
C.  Authors:
C.
C.     J. Engel
C.     T.K Gaisser
C.     P.Lipari
C.     T. Stanev
C. 
C.  This set of routines  when used in  the simulation of cosmic ray
C.  showers have only three  "contact points" with the "external world"
C.
C.    (i) SUBROUTINE NUC_NUC_INI
C.        (no  calling arguments)         
C.         to be called once during general initialization
C.    (ii) SUBROUTINE HEAVY (IA, E0)
C.         where IA (integer) is the mass number of the projectile
C.         nucleus  and E0 (TeV) is the energy per nucleon
C.         The output (positions of first interaction for the IA
C.         nucleons of the projectile) is  contained in  the common block:
C.           COMMON /C1STNC/ XX0(60),XX(60),YY(60),AX(60),AY(60)
C.         In detail:
C.             XX0(j)   (g cm-2) =  position of interaction
C.             XX(j) (mm)    x-distance from shower axis
C.             YY(j) (mm)    y-distance from shower axis
C.             AX(j) (radiants)  Theta_x with respect to original direction
C.             AY(j) (radiants)  Theta_y with respect to original direction
C.      
C.    (iii)  FUNCTION FPNI (E,L)
C.           Interaction length in air.
C.           E (TeV) is the energy of the particle, L is the particle
C.           code (NOTE: "Sibyll" codes are used : L =1-18) 
C.           WANRNING : The nucleus-nucleus cross section
C.           tabulated in the program are "matched" to the p-Air
C.           cross section calculated  with this FUNCTION, in other words 
C.           they are both calculated with the same input pp cross section
C==========================================================================

      SUBROUTINE NUC_NUC_INI

C-----------------------------------------------------------------------
C...Initialization for the generation of nucleus-nucleus interactions
C.  INPUT : E0 (TeV) Energy per nucleon of the beam nucleus
C........................................................................
      SAVE

      CALL NUC_GEOM_INI                       ! nucleus profiles
      CALL SIGMA_INI                          ! initialize pp cross sections

      RETURN
      END
C=======================================================================

      SUBROUTINE FRAGM1 (IA,NW, NF, IAF)

C-----------------------------------------------------------------------
C...Nuclear Fragmentation 
C.  total dissolution of nucleus
C.......................................................................
      SAVE

      DIMENSION IAF(60)
      NF = IA-NW
      DO J=1,NF
         IAF(J) = 1
      ENDDO
      RETURN
      END
C->
C=======================================================================

      SUBROUTINE FRAGM2 (IA,NW, NF, IAF)

C-----------------------------------------------------------------------
C...Nuclear Fragmentation 
C.  Spectator in one single fragment 
C.......................................................................
      SAVE

      DIMENSION IAF(60)
      IF (IA-NW .GT. 0)  THEN
         NF = 1
         IAF(1) = IA-NW
      ELSE
         NF = 0
      ENDIF
      RETURN
      END

C-----------------------------------------------------------------------
C...Code of fragmentation  of spectator nucleons
C.  based on Jon Engel  abrasion-ablation algorithms
C=======================================================================

      BLOCK DATA FRAG_DATA

C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)

C...Data for the fragmentation of  nucleus  projectiles
      COMMON /FRAGMOD/A(10,10,20),AE(10,10,20),ERES(10,10),NFLAGG(10,10)
      SAVE
      DATA (NFLAGG(I, 1),I=1,10)  / 
     +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
      DATA (NFLAGG(I, 2),I=1,10)  / 
     +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
      DATA (NFLAGG(I, 3),I=1,10)  / 
     +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
      DATA (NFLAGG(I, 4),I=1,10)  / 
     +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
      DATA (NFLAGG(I, 5),I=1,10)  / 
     +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
      DATA (NFLAGG(I, 6),I=1,10)  / 
     +    0,  0,  0,  0,  0,  0,  0,  1,  1,  1 /
      DATA (NFLAGG(I, 7),I=1,10)  / 
     +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
      DATA (NFLAGG(I, 8),I=1,10)  / 
     +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
      DATA (NFLAGG(I, 9),I=1,10)  / 
     +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
      DATA (NFLAGG(I,10),I=1,10)  / 
     +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
      DATA (A(I, 1, 1),I=1,10)  / 
     +  .438D-01,.172D0  ,.283D0  ,.511D0  ,.715D0  ,.920D0  ,1.19D0  ,
     +  1.37D0  ,1.65D0  ,2.14D0   /
      DATA (A(I, 1, 2),I=1,10)  / 
     +  .147D-01,.249D-01,.439D-01,.592D-01,.776D-01,.886D-01,.108D0  ,
     +  .117D0  ,.126D0  ,.128D0   /
      DATA (A(I, 1, 3),I=1,10)  / 
     +  .216D-02,.627D-02,.834D-02,.108D-01,.144D-01,.152D-01,.196D-01,
     +  .200D-01,.210D-01,.224D-01 /
      DATA (A(I, 1, 4),I=1,10)  / 
     +  .593D-01,.653D-01,.116D0  ,.145D0  ,.184D0  ,.204D0  ,.234D0  ,
     +  .257D0  ,.271D0  ,.248D0   /
      DATA (A(I, 1, 5),I=1,10)  / 
     +  .000D+00,.918D-02,.362D-02,.805D-02,.436D-02,.728D-02,.466D-02,
     +  .707D-02,.932D-02,.130D-01 /
      DATA (A(I, 1, 6),I=1,10)  / 
     +  .000D+00,.180D-02,.247D-02,.208D-02,.224D-02,.214D-02,.226D-02,
     +  .233D-02,.230D-02,.194D-02 /
      DATA (A(I, 1, 7),I=1,10)  / 
     +  .000D+00,.106D-02,.703D-03,.687D-03,.739D-03,.674D-03,.819D-03,
     +  .768D-03,.756D-03,.720D-03 /
      DATA (A(I, 1, 8),I=1,10)  / 
     +  .000D+00,.000D+00,.188D-02,.130D-02,.138D-02,.117D-02,.124D-02,
     +  .119D-02,.111D-02,.829D-03 /
      DATA (A(I, 1, 9),I=1,10)  / 
     +  .000D+00,.000D+00,.302D-03,.258D-03,.249D-03,.208D-03,.248D-03,
     +  .222D-03,.210D-03,.187D-03 /
      DATA (A(I, 1,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.235D-03,.222D-03,.172D-03,.181D-03,
     +  .166D-03,.152D-03,.124D-03 /
      DATA (A(I, 1,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.238D-03,.179D-03,.145D-03,.156D-03,
     +  .138D-03,.129D-03,.111D-03 /
      DATA (A(I, 1,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.368D-03,.400D-03,.255D-03,.262D-03,
     +  .221D-03,.182D-03,.112D-03 /
      DATA (A(I, 1,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.753D-04,.712D-04,.527D-04,
     +  .537D-04,.538D-04,.487D-04 /
      DATA (A(I, 1,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.103D-03,.589D-04,.578D-04,
     +  .468D-04,.385D-04,.269D-04 /
      DATA (A(I, 1,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.444D-04,.372D-04,
     +  .318D-04,.284D-04,.218D-04 /
      DATA (A(I, 1,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.487D-04,.473D-04,
     +  .338D-04,.243D-04,.122D-04 /
      DATA (A(I, 1,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.121D-04,.117D-04,
     +  .932D-05,.792D-05,.583D-05 /
      DATA (A(I, 1,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.147D-04,
     +  .101D-04,.756D-05,.496D-05 /
      DATA (A(I, 1,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.755D-05,
     +  .612D-05,.505D-05,.341D-05 /
      DATA (A(I, 1,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  .630D-05,.444D-05,.282D-05 /
      DATA (A(I, 2, 1),I=1,10)  / 
     +  .269D0  ,.510D0  ,.738D0  ,1.12D0  ,1.46D0  ,1.83D0  ,2.22D0  ,
     +  2.57D0  ,3.00D0  ,3.67D0   /
      DATA (A(I, 2, 2),I=1,10)  / 
     +  .121D0  ,.133D0  ,.190D0  ,.234D0  ,.293D0  ,.332D0  ,.395D0  ,
     +  .431D0  ,.468D0  ,.502D0   /
      DATA (A(I, 2, 3),I=1,10)  / 
     +  .227D-01,.374D-01,.474D-01,.578D-01,.722D-01,.794D-01,.960D-01,
     +  .102D0  ,.110D0  ,.120D0   /
      DATA (A(I, 2, 4),I=1,10)  / 
     +  .287D0  ,.196D0  ,.270D0  ,.314D0  ,.373D0  ,.408D0  ,.462D0  ,
     +  .498D0  ,.529D0  ,.523D0   /
      DATA (A(I, 2, 5),I=1,10)  / 
     +  .000D+00,.433D-01,.218D-01,.384D-01,.263D-01,.385D-01,.298D-01,
     +  .405D-01,.504D-01,.671D-01 /
      DATA (A(I, 2, 6),I=1,10)  / 
     +  .000D+00,.151D-01,.177D-01,.159D-01,.173D-01,.173D-01,.187D-01,
     +  .196D-01,.201D-01,.191D-01 /
      DATA (A(I, 2, 7),I=1,10)  / 
     +  .000D+00,.457D-02,.607D-02,.610D-02,.677D-02,.670D-02,.784D-02,
     +  .787D-02,.806D-02,.803D-02 /
      DATA (A(I, 2, 8),I=1,10)  / 
     +  .000D+00,.000D+00,.702D-02,.536D-02,.558D-02,.510D-02,.554D-02,
     +  .546D-02,.538D-02,.489D-02 /
      DATA (A(I, 2, 9),I=1,10)  / 
     +  .000D+00,.000D+00,.190D-02,.199D-02,.205D-02,.191D-02,.221D-02,
     +  .214D-02,.213D-02,.204D-02 /
      DATA (A(I, 2,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.226D-02,.219D-02,.195D-02,.208D-02,
     +  .204D-02,.203D-02,.194D-02 /
      DATA (A(I, 2,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.213D-02,.195D-02,.175D-02,.191D-02,
     +  .183D-02,.179D-02,.166D-02 /
      DATA (A(I, 2,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.588D-03,.186D-02,.137D-02,.141D-02,
     +  .128D-02,.117D-02,.947D-03 /
      DATA (A(I, 2,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.554D-03,.562D-03,.454D-03,
     +  .485D-03,.505D-03,.509D-03 /
      DATA (A(I, 2,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.490D-03,.533D-03,.531D-03,
     +  .476D-03,.437D-03,.369D-03 /
      DATA (A(I, 2,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.427D-03,.382D-03,
     +  .358D-03,.340D-03,.294D-03 /
      DATA (A(I, 2,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.239D-03,.298D-03,
     +  .238D-03,.196D-03,.134D-03 /
      DATA (A(I, 2,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.299D-04,.893D-04,
     +  .796D-04,.744D-04,.683D-04 /
      DATA (A(I, 2,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.127D-03,
     +  .107D-03,.916D-04,.720D-04 /
      DATA (A(I, 2,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.397D-04,
     +  .630D-04,.565D-04,.461D-04 /
      DATA (A(I, 2,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  .511D-04,.459D-04,.402D-04 /
      DATA (A(I, 3, 1),I=1,10)  / 
     +  .708D0  ,1.02D0  ,1.41D0  ,1.91D0  ,2.42D0  ,3.00D0  ,3.53D0  ,
     +  4.09D0  ,4.71D0  ,5.57D0   /
      DATA (A(I, 3, 2),I=1,10)  / 
     +  .397D0  ,.410D0  ,.539D0  ,.648D0  ,.795D0  ,.910D0  ,1.06D0  ,
     +  1.17D0  ,1.29D0  ,1.42D0   /
      DATA (A(I, 3, 3),I=1,10)  / 
     +  .845D-01,.122D0  ,.157D0  ,.190D0  ,.232D0  ,.262D0  ,.307D0  ,
     +  .335D0  ,.366D0  ,.402D0   /
      DATA (A(I, 3, 4),I=1,10)  / 
     +  .210D0  ,.379D0  ,.450D0  ,.490D0  ,.574D0  ,.636D0  ,.709D0  ,
     +  .769D0  ,.820D0  ,.849D0   /
      DATA (A(I, 3, 5),I=1,10)  / 
     +  .000D+00,.102D0  ,.675D-01,.104D0  ,.858D-01,.115D0  ,.102D0  ,
     +  .129D0  ,.154D0  ,.194D0   /
      DATA (A(I, 3, 6),I=1,10)  / 
     +  .000D+00,.392D-01,.615D-01,.593D-01,.649D-01,.674D-01,.735D-01,
     +  .779D-01,.817D-01,.828D-01 /
      DATA (A(I, 3, 7),I=1,10)  / 
     +  .000D+00,.539D-02,.222D-01,.238D-01,.269D-01,.280D-01,.320D-01,
     +  .334D-01,.350D-01,.361D-01 /
      DATA (A(I, 3, 8),I=1,10)  / 
     +  .000D+00,.000D+00,.838D-02,.130D-01,.133D-01,.131D-01,.141D-01,
     +  .144D-01,.149D-01,.152D-01 /
      DATA (A(I, 3, 9),I=1,10)  / 
     +  .000D+00,.000D+00,.228D-02,.647D-02,.688D-02,.687D-02,.772D-02,
     +  .786D-02,.811D-02,.824D-02 /
      DATA (A(I, 3,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.664D-02,.828D-02,.802D-02,.845D-02,
     +  .869D-02,.902D-02,.930D-02 /
      DATA (A(I, 3,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.338D-02,.735D-02,.710D-02,.767D-02,
     +  .767D-02,.776D-02,.756D-02 /
      DATA (A(I, 3,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.280D-03,.262D-02,.349D-02,.342D-02,
     +  .322D-02,.312D-02,.291D-02 /
      DATA (A(I, 3,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.618D-03,.161D-02,.138D-02,
     +  .148D-02,.155D-02,.166D-02 /
      DATA (A(I, 3,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.313D-03,.128D-02,.161D-02,
     +  .150D-02,.144D-02,.134D-02 /
      DATA (A(I, 3,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.645D-03,.118D-02,
     +  .115D-02,.111D-02,.103D-02 /
      DATA (A(I, 3,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.117D-03,.497D-03,
     +  .581D-03,.501D-03,.401D-03 /
      DATA (A(I, 3,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.115D-04,.997D-04,
     +  .202D-03,.203D-03,.206D-03 /
      DATA (A(I, 3,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.877D-04,
     +  .242D-03,.263D-03,.226D-03 /
      DATA (A(I, 3,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.158D-04,
     +  .881D-04,.152D-03,.136D-03 /
      DATA (A(I, 3,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  .358D-04,.997D-04,.117D-03 /
      DATA (A(I, 4, 1),I=1,10)  / 
     +  .945D0  ,1.29D0  ,1.40D0  ,1.98D0  ,2.73D0  ,3.17D0  ,3.77D0  ,
     +  4.29D0  ,4.78D0  ,5.54D0   /
      DATA (A(I, 4, 2),I=1,10)  / 
     +  .581D0  ,.599D0  ,.645D0  ,.839D0  ,1.10D0  ,1.25D0  ,1.47D0  ,
     +  1.64D0  ,1.78D0  ,1.99D0   /
      DATA (A(I, 4, 3),I=1,10)  / 
     +  .127D0  ,.182D0  ,.202D0  ,.264D0  ,.344D0  ,.387D0  ,.455D0  ,
     +  .504D0  ,.549D0  ,.611D0   /
      DATA (A(I, 4, 4),I=1,10)  / 
     +  .183D0  ,.464D0  ,.351D0  ,.444D0  ,.642D0  ,.659D0  ,.772D0  ,
     +  .830D0  ,.882D0  ,.930D0   /
      DATA (A(I, 4, 5),I=1,10)  / 
     +  .000D+00,.122D0  ,.803D-01,.136D0  ,.134D0  ,.173D0  ,.164D0  ,
     +  .203D0  ,.239D0  ,.300D0   /
      DATA (A(I, 4, 6),I=1,10)  / 
     +  .000D+00,.393D-01,.766D-01,.872D-01,.108D0  ,.111D0  ,.123D0  ,
     +  .132D0  ,.139D0  ,.145D0   /
      DATA (A(I, 4, 7),I=1,10)  / 
     +  .000D+00,.416D-02,.289D-01,.360D-01,.454D-01,.477D-01,.549D-01,
     +  .583D-01,.618D-01,.654D-01 /
      DATA (A(I, 4, 8),I=1,10)  / 
     +  .000D+00,.000D+00,.761D-02,.157D-01,.214D-01,.205D-01,.233D-01,
     +  .241D-01,.255D-01,.271D-01 /
      DATA (A(I, 4, 9),I=1,10)  / 
     +  .000D+00,.000D+00,.238D-02,.803D-02,.123D-01,.123D-01,.140D-01,
     +  .145D-01,.153D-01,.160D-01 /
      DATA (A(I, 4,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.695D-02,.150D-01,.154D-01,.166D-01,
     +  .172D-01,.181D-01,.192D-01 /
      DATA (A(I, 4,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.355D-02,.104D-01,.143D-01,.156D-01,
     +  .158D-01,.164D-01,.165D-01 /
      DATA (A(I, 4,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.112D-03,.276D-02,.568D-02,.736D-02,
     +  .684D-02,.691D-02,.661D-02 /
      DATA (A(I, 4,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.740D-03,.222D-02,.339D-02,
     +  .352D-02,.382D-02,.409D-02 /
      DATA (A(I, 4,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.369D-03,.160D-02,.322D-02,
     +  .375D-02,.375D-02,.355D-02 /
      DATA (A(I, 4,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.750D-03,.190D-02,
     +  .298D-02,.319D-02,.299D-02 /
      DATA (A(I, 4,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.260D-03,.673D-03,
     +  .117D-02,.156D-02,.126D-02 /
      DATA (A(I, 4,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.283D-05,.131D-03,
     +  .363D-03,.618D-03,.690D-03 /
      DATA (A(I, 4,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.205D-03,
     +  .378D-03,.709D-03,.844D-03 /
      DATA (A(I, 4,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.654D-05,
     +  .150D-03,.341D-03,.527D-03 /
      DATA (A(I, 4,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  .957D-04,.197D-03,.406D-03 /
      DATA (A(I, 5, 1),I=1,10)  / 
     +  1.16D0  ,1.70D0  ,2.19D0  ,2.79D0  ,3.33D0  ,3.90D0  ,4.49D0  ,
     +  5.07D0  ,5.66D0  ,6.38D0   /
      DATA (A(I, 5, 2),I=1,10)  / 
     +  .779D0  ,.899D0  ,1.09D0  ,1.28D0  ,1.51D0  ,1.71D0  ,1.96D0  ,
     +  2.18D0  ,2.39D0  ,2.62D0   /
      DATA (A(I, 5, 3),I=1,10)  / 
     +  .167D0  ,.263D0  ,.334D0  ,.408D0  ,.482D0  ,.548D0  ,.632D0  ,
     +  .700D0  ,.767D0  ,.840D0   /
      DATA (A(I, 5, 4),I=1,10)  / 
     +  .203D0  ,.565D0  ,.845D0  ,.867D0  ,.906D0  ,.961D0  ,1.08D0  ,
     +  1.13D0  ,1.21D0  ,1.25D0   /
      DATA (A(I, 5, 5),I=1,10)  / 
     +  .000D+00,.129D0  ,.152D0  ,.237D0  ,.208D0  ,.268D0  ,.258D0  ,
     +  .312D0  ,.368D0  ,.450D0   /
      DATA (A(I, 5, 6),I=1,10)  / 
     +  .000D+00,.460D-01,.126D0  ,.174D0  ,.182D0  ,.188D0  ,.208D0  ,
     +  .219D0  ,.233D0  ,.239D0   /
      DATA (A(I, 5, 7),I=1,10)  / 
     +  .000D+00,.289D-02,.380D-01,.611D-01,.788D-01,.845D-01,.974D-01,
     +  .103D0  ,.111D0  ,.117D0   /
      DATA (A(I, 5, 8),I=1,10)  / 
     +  .000D+00,.000D+00,.137D-01,.223D-01,.374D-01,.436D-01,.488D-01,
     +  .488D-01,.524D-01,.547D-01 /
      DATA (A(I, 5, 9),I=1,10)  / 
     +  .000D+00,.000D+00,.162D-02,.114D-01,.198D-01,.263D-01,.315D-01,
     +  .323D-01,.348D-01,.364D-01 /
      DATA (A(I, 5,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.149D-01,.240D-01,.320D-01,.428D-01,
     +  .436D-01,.469D-01,.493D-01 /
      DATA (A(I, 5,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.562D-02,.194D-01,.290D-01,.408D-01,
     +  .460D-01,.492D-01,.500D-01 /
      DATA (A(I, 5,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.476D-04,.106D-01,.134D-01,.191D-01,
     +  .227D-01,.264D-01,.253D-01 /
      DATA (A(I, 5,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.281D-02,.679D-02,.879D-02,
     +  .123D-01,.165D-01,.190D-01 /
      DATA (A(I, 5,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.542D-04,.847D-02,.125D-01,
     +  .144D-01,.173D-01,.192D-01 /
      DATA (A(I, 5,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.652D-02,.982D-02,
     +  .129D-01,.159D-01,.192D-01 /
      DATA (A(I, 5,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.109D-03,.688D-02,
     +  .751D-02,.845D-02,.905D-02 /
      DATA (A(I, 5,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.823D-06,.237D-02,
     +  .318D-02,.446D-02,.569D-02 /
      DATA (A(I, 5,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.604D-03,
     +  .610D-02,.673D-02,.827D-02 /
      DATA (A(I, 5,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.716D-06,
     +  .412D-02,.519D-02,.617D-02 /
      DATA (A(I, 5,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  .710D-03,.543D-02,.674D-02 /
      DATA (A(I, 6, 1),I=1,10)  / 
     +  1.36D0  ,2.08D0  ,2.67D0  ,3.30D0  ,3.94D0  ,4.62D0  ,5.18D0  ,
     +  3.60D0  ,3.64D0  ,3.95D0   /
      DATA (A(I, 6, 2),I=1,10)  / 
     +  1.07D0  ,1.33D0  ,1.58D0  ,1.82D0  ,2.10D0  ,2.44D0  ,2.74D0  ,
     +  1.78D0  ,1.73D0  ,1.80D0   /
      DATA (A(I, 6, 3),I=1,10)  / 
     +  .158D0  ,.276D0  ,.402D0  ,.506D0  ,.609D0  ,.700D0  ,.802D0  ,
     +  .638D0  ,.629D0  ,.658D0   /
      DATA (A(I, 6, 4),I=1,10)  / 
     +  .308D0  ,.739D0  ,1.02D0  ,1.12D0  ,1.26D0  ,1.35D0  ,1.57D0  ,
     +  1.94D0  ,1.71D0  ,1.55D0   /
      DATA (A(I, 6, 5),I=1,10)  / 
     +  .000D+00,.217D0  ,.183D0  ,.324D0  ,.276D0  ,.395D0  ,.393D0  ,
     +  .558D0  ,.602D0  ,.681D0   /
      DATA (A(I, 6, 6),I=1,10)  / 
     +  .000D+00,.658D-01,.251D0  ,.267D0  ,.299D0  ,.326D0  ,.386D0  ,
     +  .452D0  ,.475D0  ,.409D0   /
      DATA (A(I, 6, 7),I=1,10)  / 
     +  .000D+00,.198D-02,.774D-01,.136D0  ,.149D0  ,.164D0  ,.187D0  ,
     +  .210D0  ,.238D0  ,.256D0   /
      DATA (A(I, 6, 8),I=1,10)  / 
     +  .000D+00,.000D+00,.290D-01,.122D0  ,.139D0  ,.128D0  ,.129D0  ,
     +  .137D0  ,.147D0  ,.167D0   /
      DATA (A(I, 6, 9),I=1,10)  / 
     +  .000D+00,.000D+00,.699D-03,.617D-01,.750D-01,.801D-01,.905D-01,
     +  .974D-01,.105D0  ,.122D0   /
      DATA (A(I, 6,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.310D-01,.112D0  ,.127D0  ,.140D0  ,
     +  .143D0  ,.155D0  ,.176D0   /
      DATA (A(I, 6,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.277D-02,.889D-01,.143D0  ,.150D0  ,
     +  .175D0  ,.184D0  ,.208D0   /
      DATA (A(I, 6,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.202D-04,.343D-01,.959D-01,.109D0  ,
     +  .115D0  ,.112D0  ,.116D0   /
      DATA (A(I, 6,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.186D-02,.435D-01,.512D-01,
     +  .744D-01,.856D-01,.103D0   /
      DATA (A(I, 6,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.144D-04,.427D-01,.786D-01,
     +  .911D-01,.993D-01,.108D0   /
      DATA (A(I, 6,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.466D-02,.518D-01,
     +  .848D-01,.109D0  ,.119D0   /
      DATA (A(I, 6,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.655D-05,.330D-01,
     +  .586D-01,.617D-01,.594D-01 /
      DATA (A(I, 6,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.228D-06,.328D-02,
     +  .190D-01,.301D-01,.454D-01 /
      DATA (A(I, 6,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.218D-04,
     +  .272D-01,.501D-01,.707D-01 /
      DATA (A(I, 6,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.146D-06,
     +  .441D-02,.378D-01,.556D-01 /
      DATA (A(I, 6,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  .160D-03,.204D-01,.679D-01 /
      DATA (A(I, 7, 1),I=1,10)  / 
     +  .522D0  ,.862D0  ,1.14D0  ,1.40D0  ,1.70D0  ,1.94D0  ,2.26D0  ,
     +  2.48D0  ,2.72D0  ,3.95D0   /
      DATA (A(I, 7, 2),I=1,10)  / 
     +  .314D0  ,.450D0  ,.588D0  ,.692D0  ,.834D0  ,.936D0  ,1.09D0  ,
     +  1.18D0  ,1.28D0  ,1.80D0   /
      DATA (A(I, 7, 3),I=1,10)  / 
     +  .814D-01,.147D0  ,.189D0  ,.226D0  ,.272D0  ,.302D0  ,.351D0  ,
     +  .378D0  ,.406D0  ,.658D0   /
      DATA (A(I, 7, 4),I=1,10)  / 
     +  .252D0  ,.864D0  ,1.01D0  ,.851D0  ,.837D0  ,.774D0  ,.763D0  ,
     +  .757D0  ,.748D0  ,1.55D0   /
      DATA (A(I, 7, 5),I=1,10)  / 
     +  .000D+00,.225D0  ,.180D0  ,.276D0  ,.193D0  ,.240D0  ,.190D0  ,
     +  .228D0  ,.259D0  ,.681D0   /
      DATA (A(I, 7, 6),I=1,10)  / 
     +  .000D+00,.485D-01,.272D0  ,.273D0  ,.253D0  ,.216D0  ,.206D0  ,
     +  .197D0  ,.191D0  ,.409D0   /
      DATA (A(I, 7, 7),I=1,10)  / 
     +  .000D+00,.137D-02,.752D-01,.137D0  ,.152D0  ,.134D0  ,.125D0  ,
     +  .119D0  ,.116D0  ,.256D0   /
      DATA (A(I, 7, 8),I=1,10)  / 
     +  .000D+00,.000D+00,.220D-01,.155D0  ,.175D0  ,.155D0  ,.116D0  ,
     +  .977D-01,.858D-01,.167D0   /
      DATA (A(I, 7, 9),I=1,10)  / 
     +  .000D+00,.000D+00,.326D-03,.695D-01,.881D-01,.106D0  ,.897D-01,
     +  .782D-01,.706D-01,.122D0   /
      DATA (A(I, 7,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.261D-01,.124D0  ,.131D0  ,.156D0  ,
     +  .141D0  ,.121D0  ,.176D0   /
      DATA (A(I, 7,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.785D-03,.864D-01,.130D0  ,.170D0  ,
     +  .182D0  ,.172D0  ,.208     /
      DATA (A(I, 7,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.896D-05,.225D-01,.105D0  ,.126D0  ,
     +  .126D0  ,.135D0  ,.116D0   /
      DATA (A(I, 7,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.542D-03,.427D-01,.553D-01,
     +  .744D-01,.980D-01,.103D0   /
      DATA (A(I, 7,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.515D-05,.377D-01,.831D-01,
     +  .985D-01,.104D0  ,.108D0   /
      DATA (A(I, 7,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.285D-02,.495D-01,
     +  .871D-01,.106D0  ,.119D0   /
      DATA (A(I, 7,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.110D-05,.284D-01,
     +  .588D-01,.657D-01,.594D-01 /
      DATA (A(I, 7,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.722D-07,.176D-02,
     +  .170D-01,.305D-01,.454D-01 /
      DATA (A(I, 7,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.148D-05,
     +  .213D-01,.492D-01,.707D-01 /
      DATA (A(I, 7,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.323D-07,
     +  .722D-02,.359D-01,.556D-01 /
      DATA (A(I, 7,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  .461D-05,.155D-01,.679D-01 /
      DATA (A(I, 8, 1),I=1,10)  / 
     +  .630D0  ,.974D0  ,1.29D0  ,1.58D0  ,1.89D0  ,2.16D0  ,2.49D0  ,
     +  2.75D0  ,3.02D0  ,3.95D0   /
      DATA (A(I, 8, 2),I=1,10)  / 
     +  .328D0  ,.459D0  ,.613D0  ,.735D0  ,.879D0  ,.994D0  ,1.15D0  ,
     +  1.27D0  ,1.38D0  ,1.80D0   /
      DATA (A(I, 8, 3),I=1,10)  / 
     +  .748D-01,.121D0  ,.164D0  ,.197D0  ,.235D0  ,.265D0  ,.310D0  ,
     +  .339D0  ,.370D0  ,.658D0   /
      DATA (A(I, 8, 4),I=1,10)  / 
     +  .194D0  ,.211D0  ,.337D0  ,.344D0  ,.339D0  ,.351D0  ,.390    ,
     +  .419D0  ,.442D0  ,1.55D0   /
      DATA (A(I, 8, 5),I=1,10)  / 
     +  .000D+00,.869D-01,.725D-01,.113D0  ,.810D-01,.106D0  ,.951D-01,
     +  .120D0  ,.143D0  ,.681D0   /
      DATA (A(I, 8, 6),I=1,10)  / 
     +  .000D+00,.288D-01,.102D0  ,.922D-01,.857D-01,.845D-01,.932D-01,
     +  .983D-01,.102D0  ,.409D0   /
      DATA (A(I, 8, 7),I=1,10)  / 
     +  .000D+00,.668D-03,.533D-01,.575D-01,.493D-01,.482D-01,.539D-01,
     +  .558D-01,.582D-01,.256D0   /
      DATA (A(I, 8, 8),I=1,10)  / 
     +  .000D+00,.000D+00,.205D-01,.808D-01,.510D-01,.409D-01,.406D-01,
     +  .394D-01,.389D-01,.167D0   /
      DATA (A(I, 8, 9),I=1,10)  / 
     +  .000D+00,.000D+00,.999D-04,.647D-01,.385D-01,.325D-01,.325D-01,
     +  .316D-01,.314D-01,.122D0   /
      DATA (A(I, 8,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.169D-01,.834D-01,.611D-01,.565D-01,
     +  .533D-01,.519D-01,.176D0   /
      DATA (A(I, 8,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.107D-03,.769D-01,.922D-01,.805D-01,
     +  .745D-01,.711D-01,.208D0   /
      DATA (A(I, 8,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.180D-05,.143D-01,.983D-01,.775D-01,
     +  .627D-01,.541D-01,.116D0   /
      DATA (A(I, 8,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.157D-04,.346D-01,.507D-01,
     +  .479D-01,.455D-01,.103D0   /
      DATA (A(I, 8,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.752D-06,.248D-01,.721D-01,
     +  .728D-01,.611D-01,.108D0   /
      DATA (A(I, 8,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.686D-04,.356D-01,
     +  .731D-01,.791D-01,.119D0   /
      DATA (A(I, 8,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.838D-07,.151D-01,
     +  .470D-01,.567D-01,.594D-01 /
      DATA (A(I, 8,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.759D-08,.400D-04,
     +  .193D-01,.313D-01,.454D-01 /
      DATA (A(I, 8,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.385D-07,
     +  .921D-02,.353D-01,.707D-01 /
      DATA (A(I, 8,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.219D-08,
     +  .348D-03,.226D-01,.556D-01 /
      DATA (A(I, 8,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  .212D-07,.149D-01,.679D-01 /
      DATA (A(I, 9, 1),I=1,10)  / 
     +  .736D0  ,1.13D0  ,1.49D0  ,1.82D0  ,2.20D0  ,2.49D0  ,2.86D0  ,
     +  3.17D0  ,3.49D0  ,3.95D0   /
      DATA (A(I, 9, 2),I=1,10)  / 
     +  .339D0  ,.492D0  ,.658D0  ,.789D0  ,.958D0  ,1.08D0  ,1.25D0  ,
     +  1.37D0  ,1.50D0  ,1.80D0   /
      DATA (A(I, 9, 3),I=1,10)  / 
     +  .680D-01,.110D0  ,.150D0  ,.180D0  ,.222D0  ,.247D0  ,.289    ,
     +  .318D0  ,.349D0  ,.658D0   /
      DATA (A(I, 9, 4),I=1,10)  / 
     +  .110D0  ,.104D0  ,.157D0  ,.156D0  ,.210D0  ,.205D0  ,.246D0  ,
     +  .274D0  ,.300D0  ,1.55D0   /
      DATA (A(I, 9, 5),I=1,10)  / 
     +  .000D+00,.379D-01,.347D-01,.477D-01,.486D-01,.576D-01,.569D-01,
     +  .732D-01,.893D-01,.681D0   /
      DATA (A(I, 9, 6),I=1,10)  / 
     +  .000D+00,.223D-01,.354D-01,.312D-01,.436D-01,.400D-01,.489D-01,
     +  .548D-01,.600D-01,.409D0   /
      DATA (A(I, 9, 7),I=1,10)  / 
     +  .000D+00,.338D-03,.149D-01,.142D-01,.215D-01,.188D-01,.248D-01,
     +  .278D-01,.307D-01,.256D0   /
      DATA (A(I, 9, 8),I=1,10)  / 
     +  .000D+00,.000D+00,.553D-02,.862D-02,.150D-01,.106D-01,.145D-01,
     +  .165D-01,.181D-01,.167D0   /
      DATA (A(I, 9, 9),I=1,10)  / 
     +  .000D+00,.000D+00,.375D-04,.641D-02,.111D-01,.792D-02,.112D-01,
     +  .127D-01,.140D-01,.122D0   /
      DATA (A(I, 9,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.112D-01,.200D-01,.127D-01,.176D-01,
     +  .200D-01,.220D-01,.176D0   /
      DATA (A(I, 9,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.244D-04,.261D-01,.162D-01,.232D-01,
     +  .263D-01,.287D-01,.208D0   /
      DATA (A(I, 9,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.455D-06,.635D-02,.121D-01,.186D-01,
     +  .201D-01,.207D-01,.116D0   /
      DATA (A(I, 9,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.146D-05,.922D-02,.116D-01,
     +  .145D-01,.165D-01,.103D0   /
      DATA (A(I, 9,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.135D-06,.128D-01,.202D-01,
     +  .215D-01,.220D-01,.108D0   /
      DATA (A(I, 9,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.237D-05,.229D-01,
     +  .259D-01,.271D-01,.119D0   /
      DATA (A(I, 9,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.100D-07,.534D-02,
     +  .210D-01,.193D-01,.594D-01 /
      DATA (A(I, 9,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.915D-09,.847D-06,
     +  .119D-01,.125D-01,.454D-01 /
      DATA (A(I, 9,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.298D-08,
     +  .101D-01,.242D-01,.707D-01 /
      DATA (A(I, 9,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.196D-09,
     +  .243D-05,.234D-01,.556D-01 /
      DATA (A(I, 9,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  .575D-09,.364D-02,.679D-01 /
      DATA (A(I,10, 1),I=1,10)  / 
     +  .959D0  ,1.46D0  ,1.92D0  ,2.34D0  ,2.80D0  ,3.24D0  ,3.64D0  ,
     +  4.05D0  ,4.48D0  ,3.95D0   /
      DATA (A(I,10, 2),I=1,10)  / 
     +  .343D0  ,.516D0  ,.692D0  ,.836D0  ,1.01D0  ,1.16D0  ,1.31D0  ,
     +  1.46D0  ,1.61D0  ,1.80D0   /
      DATA (A(I,10, 3),I=1,10)  / 
     +  .512D-01,.837D-01,.115D0  ,.138D0  ,.169D0  ,.195D0  ,.220D0  ,
     +  .245D0  ,.270D0  ,.658D0   /
      DATA (A(I,10, 4),I=1,10)  / 
     +  .274D-01,.361D-01,.510D-01,.562D-01,.703D-01,.828D-01,.877D-01,
     +  .996D-01,.111D0  ,1.55D0   /
      DATA (A(I,10, 5),I=1,10)  / 
     +  .000D+00,.850D-02,.875D-02,.118D-01,.124D-01,.170D-01,.154D-01,
     +  .194D-01,.237D-01,.681D0   /
      DATA (A(I,10, 6),I=1,10)  / 
     +  .000D+00,.345D-02,.519D-02,.533D-02,.691D-02,.842D-02,.844D-02,
     +  .987D-02,.113D-01,.409D0   /
      DATA (A(I,10, 7),I=1,10)  / 
     +  .000D+00,.722D-04,.130D-02,.135D-02,.189D-02,.240D-02,.235D-02,
     +  .281D-02,.331D-02,.256D0   /
      DATA (A(I,10, 8),I=1,10)  / 
     +  .000D+00,.000D+00,.283D-03,.272D-03,.394D-03,.557D-03,.480D-03,
     +  .616D-03,.775D-03,.167D0   /
      DATA (A(I,10, 9),I=1,10)  / 
     +  .000D+00,.000D+00,.457D-05,.122D-03,.192D-03,.275D-03,.225D-03,
     +  .292D-03,.373D-03,.122D0   /
      DATA (A(I,10,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.119D-03,.185D-03,.278D-03,.201D-03,
     +  .274D-03,.364D-03,.176D0   /
      DATA (A(I,10,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.140D-05,.129D-03,.200D-03,.137D-03,
     +  .188D-03,.252D-03,.208D0   /
      DATA (A(I,10,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.207D-07,.307D-04,.518D-04,.278D-04,
     +  .421D-04,.608D-04,.116D0   /
      DATA (A(I,10,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.306D-07,.252D-04,.111D-04,
     +  .188D-04,.295D-04,.103D0   /
      DATA (A(I,10,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.321D-08,.220D-04,.104D-04,
     +  .162D-04,.243D-04,.108D0   /
      DATA (A(I,10,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.770D-08,.632D-05,
     +  .105D-04,.162D-04,.119D0   /
      DATA (A(I,10,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.117D-09,.199D-05,
     +  .321D-05,.492D-05,.594D-01 /
      DATA (A(I,10,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.888E-11,.323D-09,
     +  .106D-05,.192D-05,.454D-01 /
      DATA (A(I,10,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.174E-10,
     +  .131D-05,.218D-05,.707D-01 /
      DATA (A(I,10,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.994E-12,
     +  .233D-09,.104D-05,.556D-01 /
      DATA (A(I,10,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  .144E-11,.724D-06,.679D-01 /
      DATA (AE(I, 1, 1),I=1,10)  / 
     +  7.27D0  ,6.29D0  ,7.76D0  ,6.70D0  ,8.17D0  ,7.34D0  ,8.70D0  ,
     +  8.02D0  ,7.37D0  ,6.18D0   /
      DATA (AE(I, 1, 2),I=1,10)  / 
     +  7.41D0  ,7.52D0  ,8.14D0  ,8.20D0  ,8.96D0  ,9.05D0  ,9.96D0  ,
     +  10.0D0  ,10.1D0  ,9.86D0   /
      DATA (AE(I, 1, 3),I=1,10)  / 
     +  7.72D0  ,7.69D0  ,9.17D0  ,8.99D0  ,10.6D0  ,10.5D0  ,12.1D0  ,
     +  12.1D0  ,12.0D0  ,11.5D0   /
      DATA (AE(I, 1, 4),I=1,10)  / 
     +  7.90D0  ,8.48D0  ,9.50D0  ,9.94D0  ,10.8D0  ,11.4D0  ,12.2D0  ,
     +  12.8D0  ,13.3D0  ,13.8D0   /
      DATA (AE(I, 1, 5),I=1,10)  / 
     +  .000D+00,8.52D0  ,9.59D0  ,10.1D0  ,11.1D0  ,11.8D0  ,12.7D0  ,
     +  13.3D0  ,13.8D0  ,14.4D0   /
      DATA (AE(I, 1, 6),I=1,10)  / 
     +  .000D+00,9.00D0  ,10.7D0  ,11.7D0  ,13.2D0  ,14.2D0  ,15.6D0  ,
     +  16.5D0  ,17.3D0  ,18.0D0   /
      DATA (AE(I, 1, 7),I=1,10)  / 
     +  .000D+00,9.01D0  ,11.1D0  ,11.9D0  ,14.3D0  ,15.0D0  ,17.4D0  ,
     +  18.0D0  ,18.6D0  ,18.8D0   /
      DATA (AE(I, 1, 8),I=1,10)  / 
     +  .000D+00,.000D+00,11.2D0  ,12.4D0  ,14.5D0  ,15.7D0  ,17.6D0  ,
     +  18.8D0  ,19.9D0  ,20.9D0   /
      DATA (AE(I, 1, 9),I=1,10)  / 
     +  .000D+00,.000D+00,11.4D0  ,12.7D0  ,15.5D0  ,16.6D0  ,19.3D0  ,
     +  20.2D0  ,21.1D0  ,21.7D0   /
      DATA (AE(I, 1,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,13.2D0  ,15.8D0  ,17.3D0  ,19.9D0  ,
     +  21.2D0  ,22.4D0  ,23.2D0   /
      DATA (AE(I, 1,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,13.2D0  ,16.3D0  ,17.8D0  ,20.8D0  ,
     +  22.1D0  ,23.3D0  ,24.2D0   /
      DATA (AE(I, 1,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,13.4D0  ,16.2D0  ,18.2D0  ,21.0D0  ,
     +  22.8D0  ,24.4D0  ,25.9D0   /
      DATA (AE(I, 1,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,16.5D0  ,18.4D0  ,21.6D0  ,
     +  23.2D0  ,24.8D0  ,26.2D0   /
      DATA (AE(I, 1,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,16.7D0  ,19.0D0  ,22.3D0  ,
     +  24.3D0  ,26.1D0  ,27.4D0   /
      DATA (AE(I, 1,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,19.1D0  ,22.8D0  ,
     +  24.7D0  ,26.6D0  ,28.2D0   /
      DATA (AE(I, 1,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,19.2D0  ,23.0D0  ,
     +  25.3D0  ,27.5D0  ,29.5D0   /
      DATA (AE(I, 1,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,19.6D0  ,23.3D0  ,
     +  25.6D0  ,27.8D0  ,29.6D0   /
      DATA (AE(I, 1,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,23.6D0  ,
     +  26.2D0  ,28.5D0  ,30.4D0   /
      DATA (AE(I, 1,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,23.7D0  ,
     +  26.3D0  ,28.8D0  ,31.0D0   /
      DATA (AE(I, 1,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  26.5D0  ,29.2D0  ,31.5D0   /
      DATA (AE(I, 2, 1),I=1,10)  / 
     +  8.74D0  ,8.16D0  ,9.25D0  ,8.45D0  ,9.46D0  ,8.90D0  ,9.83D0  ,
     +  9.38D0  ,8.96D0  ,8.15D0   /
      DATA (AE(I, 2, 2),I=1,10)  / 
     +  8.96D0  ,9.30D0  ,9.95D0  ,10.0D0  ,10.8D0  ,10.9D0  ,11.7D0  ,
     +  11.8D0  ,11.9D0  ,11.8D0   /
      DATA (AE(I, 2, 3),I=1,10)  / 
     +  9.44D0  ,9.66D0  ,11.0D0  ,11.0D0  ,12.3D0  ,12.5D0  ,13.7D0  ,
     +  13.9D0  ,14.0D0  ,13.8D0   /
      DATA (AE(I, 2, 4),I=1,10)  / 
     +  8.86D0  ,9.81D0  ,10.8D0  ,11.2D0  ,12.0D0  ,12.6D0  ,13.4D0  ,
     +  14.0D0  ,14.5D0  ,15.1D0   /
      DATA (AE(I, 2, 5),I=1,10)  / 
     +  .000D+00,10.2D0  ,11.4D0  ,12.0D0  ,12.9D0  ,13.6D0  ,14.5D0  ,
     +  15.1D0  ,15.7D0  ,16.3D0   /
      DATA (AE(I, 2, 6),I=1,10)  / 
     +  .000D+00,10.7D0  ,12.5D0  ,13.5D0  ,15.1D0  ,16.0D0  ,17.5D0  ,
     +  18.3D0  ,19.2D0  ,19.9D0   /
      DATA (AE(I, 2, 7),I=1,10)  / 
     +  .000D+00,11.5D0  ,12.9D0  ,13.9D0  ,16.1D0  ,17.0D0  ,19.1D0  ,
     +  19.8D0  ,20.6D0  ,21.0D0   /
      DATA (AE(I, 2, 8),I=1,10)  / 
     +  .000D+00,.000D+00,12.4D0  ,13.8D0  ,15.9D0  ,17.2D0  ,19.1D0  ,
     +  20.3D0  ,21.4D0  ,22.3D0   /
      DATA (AE(I, 2, 9),I=1,10)  / 
     +  .000D+00,.000D+00,13.4D0  ,14.5D0  ,17.1D0  ,18.3D0  ,20.9D0  ,
     +  21.9D0  ,23.0D0  ,23.7D0   /
      DATA (AE(I, 2,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,14.9D0  ,17.5D0  ,19.1D0  ,21.6D0  ,
     +  22.9D0  ,24.1D0  ,25.0D0   /
      DATA (AE(I, 2,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,15.0D0  ,18.0D0  ,19.6D0  ,22.4D0  ,
     +  23.8D0  ,25.2D0  ,26.2D0   /
      DATA (AE(I, 2,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,16.2D0  ,17.3D0  ,19.4D0  ,22.2D0  ,
     +  24.0D0  ,25.7D0  ,27.2D0   /
      DATA (AE(I, 2,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,17.8D0  ,19.8D0  ,22.9D0  ,
     +  24.6D0  ,26.2D0  ,27.7D0   /
      DATA (AE(I, 2,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,19.1D0  ,20.4D0  ,23.7D0  ,
     +  25.7D0  ,27.6D0  ,29.1D0   /
      DATA (AE(I, 2,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,20.5D0  ,24.1D0  ,
     +  26.1D0  ,28.1D0  ,29.9D0   /
      DATA (AE(I, 2,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,20.9D0  ,23.9D0  ,
     +  26.4D0  ,28.7D0  ,30.7D0   /
      DATA (AE(I, 2,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,22.4D0  ,24.2D0  ,
     +  26.7D0  ,29.0D0  ,30.9D0   /
      DATA (AE(I, 2,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,24.8D0  ,
     +  27.3D0  ,29.7D0  ,31.8D0   /
      DATA (AE(I, 2,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,26.1D0  ,
     +  27.3D0  ,29.9D0  ,32.3D0   /
      DATA (AE(I, 2,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  27.4D0  ,30.1D0  ,32.6D0   /
      DATA (AE(I, 3, 1),I=1,10)  / 
     +  11.0D0  ,11.0D0  ,11.7D0  ,11.3D0  ,11.9D0  ,11.4D0  ,12.1D0  ,
     +  11.7D0  ,11.5D0  ,11.0D0   /
      DATA (AE(I, 3, 2),I=1,10)  / 
     +  11.2D0  ,12.0D0  ,12.7D0  ,12.9D0  ,13.6D0  ,13.7D0  ,14.4D0  ,
     +  14.6D0  ,14.7D0  ,14.6D0   /
      DATA (AE(I, 3, 3),I=1,10)  / 
     +  12.1D0  ,12.6D0  ,13.7D0  ,13.9D0  ,15.0D0  ,15.2D0  ,16.3D0  ,
     +  16.5D0  ,16.7D0  ,16.7D0   /
      DATA (AE(I, 3, 4),I=1,10)  / 
     +  12.6D0  ,11.3D0  ,12.4D0  ,13.0D0  ,13.8D0  ,14.2D0  ,15.0D0  ,
     +  15.6D0  ,16.1D0  ,16.6D0   /
      DATA (AE(I, 3, 5),I=1,10)  / 
     +  .000D+00,12.6D0  ,13.7D0  ,14.4D0  ,15.3D0  ,16.0D0  ,16.8D0  ,
     +  17.5D0  ,18.1D0  ,18.6D0   /
      DATA (AE(I, 3, 6),I=1,10)  / 
     +  .000D+00,14.0D0  ,14.6D0  ,15.8D0  ,17.4D0  ,18.4D0  ,19.8D0  ,
     +  20.6D0  ,21.5D0  ,22.2D0   /
      DATA (AE(I, 3, 7),I=1,10)  / 
     +  .000D+00,16.0D0  ,15.2D0  ,16.3D0  ,18.3D0  ,19.3D0  ,21.1D0  ,
     +  22.0D0  ,22.8D0  ,23.5D0   /
      DATA (AE(I, 3, 8),I=1,10)  / 
     +  .000D+00,.000D+00,15.6D0  ,15.1D0  ,17.2D0  ,18.6D0  ,20.6D0  ,
     +  21.8D0  ,22.9D0  ,23.8D0   /
      DATA (AE(I, 3, 9),I=1,10)  / 
     +  .000D+00,.000D+00,17.8D0  ,16.3D0  ,18.8D0  ,20.1D0  ,22.5D0  ,
     +  23.6D0  ,24.7D0  ,25.6D0   /
      DATA (AE(I, 3,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,17.5D0  ,19.0D0  ,20.7D0  ,23.1D0  ,
     +  24.5D0  ,25.8D0  ,26.8D0   /
      DATA (AE(I, 3,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,19.2D0  ,19.4D0  ,21.1D0  ,23.8D0  ,
     +  25.4D0  ,26.8D0  ,28.0D0   /
      DATA (AE(I, 3,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,20.7D0  ,19.6D0  ,19.7D0  ,22.4D0  ,
     +  24.4D0  ,26.2D0  ,27.9D0   /
      DATA (AE(I, 3,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,21.6D0  ,20.4D0  ,23.2D0  ,
     +  25.1D0  ,26.9D0  ,28.5D0   /
      DATA (AE(I, 3,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,23.5D0  ,22.0D0  ,23.8D0  ,
     +  26.1D0  ,28.1D0  ,29.9D0   /
      DATA (AE(I, 3,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,23.7D0  ,24.2D0  ,
     +  26.3D0  ,28.5D0  ,30.4D0   /
      DATA (AE(I, 3,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,25.4D0  ,24.8D0  ,
     +  25.6D0  ,28.1D0  ,30.5D0   /
      DATA (AE(I, 3,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,26.9D0  ,26.8D0  ,
     +  26.1D0  ,28.4D0  ,30.8D0   /
      DATA (AE(I, 3,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,28.8D0  ,
     +  27.6D0  ,29.0D0  ,31.5D0   /
      DATA (AE(I, 3,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,30.5D0  ,
     +  29.2D0  ,28.9D0  ,31.5D0   /
      DATA (AE(I, 3,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  31.0D0  ,30.0D0  ,31.7D0   /
      DATA (AE(I, 4, 1),I=1,10)  / 
     +  13.0D0  ,13.2D0  ,14.8D0  ,14.2D0  ,14.2D0  ,14.1D0  ,14.5D0  ,
     +  14.4D0  ,14.3D0  ,14.0D0   /
      DATA (AE(I, 4, 2),I=1,10)  / 
     +  13.5D0  ,14.5D0  ,16.1D0  ,15.9D0  ,16.0D0  ,16.3D0  ,16.8D0  ,
     +  17.0D0  ,17.1D0  ,17.2D0   /
      DATA (AE(I, 4, 3),I=1,10)  / 
     +  14.9D0  ,15.3D0  ,17.2D0  ,17.1D0  ,17.5D0  ,17.8D0  ,18.6D0  ,
     +  18.9D0  ,19.1D0  ,19.3D0   /
      DATA (AE(I, 4, 4),I=1,10)  / 
     +  15.1D0  ,13.5D0  ,16.4D0  ,16.7D0  ,16.4D0  ,17.3D0  ,17.8D0  ,
     +  18.5D0  ,19.0D0  ,19.6D0   /
      DATA (AE(I, 4, 5),I=1,10)  / 
     +  .000D+00,15.6D0  ,17.5D0  ,17.7D0  ,17.8D0  ,18.6D0  ,19.2D0  ,
     +  19.9D0  ,20.3D0  ,21.1D0   /
      DATA (AE(I, 4, 6),I=1,10)  / 
     +  .000D+00,18.0D0  ,18.4D0  ,19.2D0  ,19.8D0  ,20.9D0  ,22.0D0  ,
     +  23.1D0  ,23.6D0  ,24.7D0   /
      DATA (AE(I, 4, 7),I=1,10)  / 
     +  .000D+00,27.4D0  ,19.1D0  ,19.8D0  ,20.7D0  ,21.8D0  ,23.2D0  ,
     +  24.4D0  ,24.9D0  ,25.9D0   /
      DATA (AE(I, 4, 8),I=1,10)  / 
     +  .000D+00,.000D+00,18.9D0  ,18.9D0  ,19.3D0  ,21.1D0  ,22.5D0  ,
     +  24.0D0  ,24.7D0  ,26.0D0   /
      DATA (AE(I, 4, 9),I=1,10)  / 
     +  .000D+00,.000D+00,21.1D0  ,19.7D0  ,20.7D0  ,22.3D0  ,24.0D0  ,
     +  25.6D0  ,26.3D0  ,27.7D0   /
      DATA (AE(I, 4,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,21.0D0  ,21.1D0  ,22.9D0  ,24.6D0  ,
     +  26.5D0  ,27.3D0  ,29.0D0   /
      DATA (AE(I, 4,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,21.3D0  ,22.4D0  ,23.1D0  ,25.0D0  ,
     +  27.1D0  ,27.9D0  ,29.8D0   /
      DATA (AE(I, 4,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,36.6D0  ,21.5D0  ,22.2D0  ,23.1D0  ,
     +  25.6D0  ,26.8D0  ,29.1D0   /
      DATA (AE(I, 4,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,22.9D0  ,23.1D0  ,23.7D0  ,
     +  26.2D0  ,27.3D0  ,29.6D0   /
      DATA (AE(I, 4,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,30.5D0  ,23.6D0  ,25.0D0  ,
     +  26.9D0  ,28.2D0  ,30.7D0   /
      DATA (AE(I, 4,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,25.4D0  ,26.2D0  ,
     +  27.2D0  ,28.3D0  ,31.0D0   /
      DATA (AE(I, 4,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,24.5D0  ,25.9D0  ,
     +  27.4D0  ,27.6D0  ,30.7D0   /
      DATA (AE(I, 4,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,43.3D0  ,28.4D0  ,
     +  27.5D0  ,27.9D0  ,30.9D0   /
      DATA (AE(I, 4,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,27.2D0  ,
     +  29.1D0  ,29.0D0  ,31.4D0   /
      DATA (AE(I, 4,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,51.3D0  ,
     +  30.6D0  ,29.5D0  ,31.4D0   /
      DATA (AE(I, 4,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  28.8D0  ,30.6D0  ,32.4D0   /
      DATA (AE(I, 5, 1),I=1,10)  / 
     +  15.0D0  ,14.9D0  ,15.5D0  ,15.4D0  ,15.9D0  ,15.8D0  ,16.2D0  ,
     +  16.2D0  ,16.1D0  ,15.9D0   /
      DATA (AE(I, 5, 2),I=1,10)  / 
     +  15.4D0  ,16.1D0  ,17.0D0  ,17.4D0  ,18.0D0  ,18.2D0  ,18.7D0  ,
     +  18.9D0  ,19.0D0  ,19.1D0   /
      DATA (AE(I, 5, 3),I=1,10)  / 
     +  17.1D0  ,17.2D0  ,18.3D0  ,18.7D0  ,19.3D0  ,19.6D0  ,20.3D0  ,
     +  20.6D0  ,20.8D0  ,20.9D0   /
      DATA (AE(I, 5, 4),I=1,10)  / 
     +  14.7D0  ,14.8D0  ,15.0D0  ,16.0D0  ,17.0D0  ,17.7D0  ,18.1D0  ,
     +  19.0D0  ,19.4D0  ,20.0D0   /
      DATA (AE(I, 5, 5),I=1,10)  / 
     +  .000D+00,16.7D0  ,17.6D0  ,18.1D0  ,18.6D0  ,19.2D0  ,19.7D0  ,
     +  20.4D0  ,20.8D0  ,21.2D0   /
      DATA (AE(I, 5, 6),I=1,10)  / 
     +  .000D+00,17.8D0  ,18.2D0  ,19.2D0  ,20.0D0  ,21.0D0  ,21.9D0  ,
     +  23.0D0  ,23.6D0  ,24.3D0   /
      DATA (AE(I, 5, 7),I=1,10)  / 
     +  .000D+00,35.2D0  ,18.9D0  ,20.3D0  ,20.6D0  ,21.5D0  ,22.6D0  ,
     +  23.7D0  ,24.2D0  ,24.7D0   /
      DATA (AE(I, 5, 8),I=1,10)  / 
     +  .000D+00,.000D+00,16.4D0  ,18.9D0  ,18.8D0  ,19.6D0  ,20.7D0  ,
     +  22.3D0  ,23.1D0  ,23.9D0   /
      DATA (AE(I, 5, 9),I=1,10)  / 
     +  .000D+00,.000D+00,33.9D0  ,19.8D0  ,20.3D0  ,20.7D0  ,21.9D0  ,
     +  23.4D0  ,24.1D0  ,24.8D0   /
      DATA (AE(I, 5,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,18.0D0  ,20.0D0  ,21.4D0  ,22.0D0  ,
     +  23.8D0  ,24.6D0  ,25.4D0   /
      DATA (AE(I, 5,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,26.4D0  ,20.4D0  ,21.2D0  ,22.3D0  ,
     +  23.8D0  ,24.7D0  ,25.5D0   /
      DATA (AE(I, 5,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,41.7D0  ,18.2D0  ,19.8D0  ,21.1D0  ,
     +  22.6D0  ,23.4D0  ,24.6D0   /
      DATA (AE(I, 5,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,22.5D0  ,20.0D0  ,21.7D0  ,
     +  22.8D0  ,23.7D0  ,24.7D0   /
      DATA (AE(I, 5,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,54.1D0  ,19.9D0  ,21.9D0  ,
     +  23.2D0  ,24.3D0  ,25.3D0   /
      DATA (AE(I, 5,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,21.2D0  ,22.2D0  ,
     +  23.6D0  ,24.9D0  ,25.5D0   /
      DATA (AE(I, 5,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,44.9D0  ,21.9D0  ,
     +  23.8D0  ,25.2D0  ,25.6D0   /
      DATA (AE(I, 5,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,47.8D0  ,22.7D0  ,
     +  23.8D0  ,24.9D0  ,26.3D0   /
      DATA (AE(I, 5,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,35.5D0  ,
     +  23.9D0  ,25.9D0  ,26.6D0   /
      DATA (AE(I, 5,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,64.3D0  ,
     +  24.1D0  ,25.7D0  ,27.1D0   /
      DATA (AE(I, 5,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  34.0D0  ,25.7D0  ,27.7D0   /
      DATA (AE(I, 6, 1),I=1,10)  / 
     +  16.6D0  ,16.5D0  ,16.8D0  ,16.7D0  ,17.0D0  ,16.5D0  ,16.7D0  ,
     +  18.3D0  ,18.9D0  ,19.0D0   /
      DATA (AE(I, 6, 2),I=1,10)  / 
     +  16.2D0  ,16.6D0  ,17.2D0  ,17.4D0  ,17.9D0  ,17.4D0  ,17.7D0  ,
     +  20.7D0  ,22.0D0  ,22.6D0   /
      DATA (AE(I, 6, 3),I=1,10)  / 
     +  18.9D0  ,18.7D0  ,18.8D0  ,18.6D0  ,18.9D0  ,18.6D0  ,18.9D0  ,
     +  21.0D0  ,22.3D0  ,22.9D0   /
      DATA (AE(I, 6, 4),I=1,10)  / 
     +  18.3D0  ,12.7D0  ,14.2D0  ,15.0D0  ,15.7D0  ,16.1D0  ,16.3D0  ,
     +  16.5D0  ,17.9D0  ,19.0D0   /
      DATA (AE(I, 6, 5),I=1,10)  / 
     +  .000D+00,15.7D0  ,15.1D0  ,15.3D0  ,16.5D0  ,16.4D0  ,16.4D0  ,
     +  17.0D0  ,18.3D0  ,19.4D0   /
      DATA (AE(I, 6, 6),I=1,10)  / 
     +  .000D+00,22.9D0  ,14.9D0  ,15.2D0  ,16.2D0  ,16.9D0  ,17.4D0  ,
     +  18.2D0  ,19.5D0  ,21.1D0   /
      DATA (AE(I, 6, 7),I=1,10)  / 
     +  .000D+00,40.7D0  ,18.4D0  ,15.9D0  ,17.1D0  ,17.7D0  ,18.9D0  ,
     +  19.5D0  ,20.3D0  ,21.1D0   /
      DATA (AE(I, 6, 8),I=1,10)  / 
     +  .000D+00,.000D+00,23.3D0  ,16.2D0  ,16.3D0  ,17.3D0  ,18.7D0  ,
     +  19.5D0  ,20.3D0  ,21.1D0   /
      DATA (AE(I, 6, 9),I=1,10)  / 
     +  .000D+00,.000D+00,49.2D0  ,19.0D0  ,19.1D0  ,19.4D0  ,20.2D0  ,
     +  20.8D0  ,21.6D0  ,22.0D0   /
      DATA (AE(I, 6,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,27.2D0  ,21.2D0  ,20.8D0  ,21.4D0  ,
     +  22.3D0  ,22.8D0  ,23.3D0   /
      DATA (AE(I, 6,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,45.6D0  ,25.0D0  ,22.8D0  ,23.9D0  ,
     +  23.6D0  ,24.3D0  ,24.4D0   /
      DATA (AE(I, 6,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,45.8D0  ,29.7D0  ,25.1D0  ,25.3D0  ,
     +  25.3D0  ,26.0D0  ,26.3D0   /
      DATA (AE(I, 6,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,42.7D0  ,29.0D0  ,28.0D0  ,
     +  27.0D0  ,27.2D0  ,27.6D0   /
      DATA (AE(I, 6,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,62.0D0  ,32.0D0  ,30.0D0  ,
     +  29.8D0  ,29.5D0  ,29.6D0   /
      DATA (AE(I, 6,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,44.5D0  ,34.4D0  ,
     +  32.7D0  ,31.5D0  ,31.8D0   /
      DATA (AE(I, 6,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,75.6D0  ,37.1D0  ,
     +  34.6D0  ,34.4D0  ,34.4D0   /
      DATA (AE(I, 6,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,51.2D0  ,45.2D0  ,
     +  39.0D0  ,37.5D0  ,36.4D0   /
      DATA (AE(I, 6,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,74.9D0  ,
     +  42.3D0  ,39.9D0  ,38.3D0   /
      DATA (AE(I, 6,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,69.5D0  ,
     +  50.7D0  ,42.3D0  ,41.4D0   /
      DATA (AE(I, 6,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  66.3D0  ,48.0D0  ,43.4D0   /
      DATA (AE(I, 7, 1),I=1,10)  / 
     +  27.0D0  ,25.8D0  ,26.3D0  ,26.2D0  ,26.7D0  ,26.7D0  ,27.1D0  ,
     +  27.1D0  ,27.2D0  ,19.0D0   /
      DATA (AE(I, 7, 2),I=1,10)  / 
     +  29.1D0  ,28.9D0  ,29.7D0  ,30.3D0  ,31.0D0  ,31.4D0  ,32.0D0  ,
     +  32.3D0  ,32.7D0  ,22.6D0   /
      DATA (AE(I, 7, 3),I=1,10)  / 
     +  31.6D0  ,29.7D0  ,30.9D0  ,31.4D0  ,32.5D0  ,33.1D0  ,34.0D0  ,
     +  34.6D0  ,35.1D0  ,22.9D0   /
      DATA (AE(I, 7, 4),I=1,10)  / 
     +  27.4D0  ,19.9D0  ,20.8D0  ,22.8D0  ,24.6D0  ,26.4D0  ,28.2D0  ,
     +  29.6D0  ,30.8D0  ,19.0D0   /
      DATA (AE(I, 7, 5),I=1,10)  / 
     +  .000D+00,24.6D0  ,24.1D0  ,25.0D0  ,27.2D0  ,28.7D0  ,30.7D0  ,
     +  31.8D0  ,32.9D0  ,19.4D0   /
      DATA (AE(I, 7, 6),I=1,10)  / 
     +  .000D+00,35.6D0  ,25.2D0  ,25.6D0  ,27.9D0  ,30.4D0  ,32.7D0  ,
     +  34.6D0  ,36.3D0  ,21.1D0   /
      DATA (AE(I, 7, 7),I=1,10)  / 
     +  .000D+00,45.4D0  ,30.9D0  ,28.2D0  ,29.0D0  ,31.2D0  ,34.0D0  ,
     +  35.8D0  ,37.4D0  ,21.1D0   /
      DATA (AE(I, 7, 8),I=1,10)  / 
     +  .000D+00,.000D+00,38.2D0  ,29.6D0  ,29.4D0  ,30.3D0  ,33.2D0  ,
     +  35.5D0  ,37.6D0  ,21.1D0   /
      DATA (AE(I, 7, 9),I=1,10)  / 
     +  .000D+00,.000D+00,59.3D0  ,34.5D0  ,33.7D0  ,32.9D0  ,35.4D0  ,
     +  37.6D0  ,39.6D0  ,22.0D0   /
      DATA (AE(I, 7,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,44.5D0  ,37.8D0  ,37.5D0  ,37.2D0  ,
     +  39.0D0  ,41.4D0  ,23.3D0   /
      DATA (AE(I, 7,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,67.0D0  ,43.6D0  ,42.0D0  ,40.8D0  ,
     +  41.4D0  ,43.0D0  ,24.4D0   /
      DATA (AE(I, 7,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,49.9D0  ,50.9D0  ,44.6D0  ,43.9D0  ,
     +  44.2D0  ,44.2D0  ,26.3D0   /
      DATA (AE(I, 7,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,67.2D0  ,50.5D0  ,48.7D0  ,
     +  48.1D0  ,47.2D0  ,27.6D0   /
      DATA (AE(I, 7,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,68.1D0  ,55.2D0  ,52.3D0  ,
     +  51.5D0  ,51.6D0  ,29.6D0   /
      DATA (AE(I, 7,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,68.7D0  ,58.6D0  ,
     +  56.5D0  ,55.7D0  ,31.8D0   /
      DATA (AE(I, 7,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,89.3D0  ,62.9D0  ,
     +  60.0D0  ,59.1D0  ,34.4D0   /
      DATA (AE(I, 7,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,56.0D0  ,72.9D0  ,
     +  66.3D0  ,64.2D0  ,36.4D0   /
      DATA (AE(I, 7,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,105.D0  ,
     +  71.3D0  ,68.3D0  ,38.3D0   /
      DATA (AE(I, 7,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,73.4D0  ,
     +  76.8D0  ,72.4D0  ,41.4D0   /
      DATA (AE(I, 7,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  107.D0  ,79.9D0  ,43.4D0   /
      DATA (AE(I, 8, 1),I=1,10)  / 
     +  35.5D0  ,35.3D0  ,35.7D0  ,35.7D0  ,36.3D0  ,36.3D0  ,36.7D0  ,
     +  36.7D0  ,36.7D0  ,19.0D0   /
      DATA (AE(I, 8, 2),I=1,10)  / 
     +  40.6D0  ,41.4D0  ,41.9D0  ,42.3D0  ,43.2D0  ,43.5D0  ,44.0D0  ,
     +  44.3D0  ,44.5D0  ,22.6D0   /
      DATA (AE(I, 8, 3),I=1,10)  / 
     +  45.4D0  ,45.7D0  ,46.4D0  ,47.0D0  ,48.1D0  ,48.7D0  ,49.4D0  ,
     +  49.8D0  ,50.2D0  ,22.9D0   /
      DATA (AE(I, 8, 4),I=1,10)  / 
     +  43.9D0  ,44.3D0  ,43.4D0  ,45.1D0  ,47.3D0  ,48.7D0  ,49.6D0  ,
     +  50.5D0  ,51.3D0  ,19.0D0   /
      DATA (AE(I, 8, 5),I=1,10)  / 
     +  .000D+00,49.3D0  ,49.6D0  ,50.5D0  ,53.2D0  ,54.2D0  ,55.4D0  ,
     +  56.1D0  ,56.8D0  ,19.4D0   /
      DATA (AE(I, 8, 6),I=1,10)  / 
     +  .000D+00,59.1D0  ,53.0D0  ,55.4D0  ,58.0D0  ,60.0D0  ,61.2D0  ,
     +  62.5D0  ,63.6D0  ,21.1D0   /
      DATA (AE(I, 8, 7),I=1,10)  / 
     +  .000D+00,54.5D0  ,57.1D0  ,59.2D0  ,62.3D0  ,64.4D0  ,66.0D0  ,
     +  67.3D0  ,68.5D0  ,21.1D0   /
      DATA (AE(I, 8, 8),I=1,10)  / 
     +  .000D+00,.000D+00,65.9D0  ,62.1D0  ,65.1D0  ,67.6D0  ,69.4D0  ,
     +  71.1D0  ,72.6D0  ,21.1D0   /
      DATA (AE(I, 8, 9),I=1,10)  / 
     +  .000D+00,.000D+00,72.2D0  ,67.1D0  ,70.5D0  ,73.1D0  ,75.1D0  ,
     +  76.8D0  ,78.4D0  ,22.0D0   /
      DATA (AE(I, 8,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,80.1D0  ,75.0D0  ,78.0D0  ,80.0D0  ,
     +  82.1D0  ,83.9D0  ,23.3D0   /
      DATA (AE(I, 8,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,94.5D0  ,82.2D0  ,82.8D0  ,85.1D0  ,
     +  87.3D0  ,89.2D0  ,24.4D0   /
      DATA (AE(I, 8,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,56.8D0  ,92.5D0  ,87.2D0  ,89.4D0  ,
     +  91.9D0  ,94.1D0  ,26.3D0   /
      DATA (AE(I, 8,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,116.D0  ,96.2D0  ,94.4D0  ,
     +  97.0D0  ,99.2D0  ,27.6D0   /
      DATA (AE(I, 8,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,78.1D0  ,104.D0  ,102.D0  ,
     +  102.D0  ,105.D0  ,29.6D0   /
      DATA (AE(I, 8,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,128.D0  ,111.D0  ,
     +  109.D0  ,110.D0  ,31.8D0   /
      DATA (AE(I, 8,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,104.D0  ,118.D0  ,
     +  117.D0  ,115.D0  ,34.4D0   /
      DATA (AE(I, 8,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,64.4D0  ,138.D0  ,
     +  124.D0  ,122.D0  ,36.4D0   /
      DATA (AE(I, 8,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,133.D0  ,
     +  133.D0  ,132.D0  ,38.3D0   /
      DATA (AE(I, 8,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,83.6D0  ,
     +  146.D0  ,139.D0  ,41.4D0   /
      DATA (AE(I, 8,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  166.D0  ,147.D0  ,43.4D0   /
      DATA (AE(I, 9, 1),I=1,10)  / 
     +  43.3D0  ,43.2D0  ,43.6D0  ,43.8D0  ,44.1D0  ,44.3D0  ,44.7D0  ,
     +  44.8D0  ,44.8D0  ,19.0D0   /
      DATA (AE(I, 9, 2),I=1,10)  / 
     +  50.9D0  ,51.4D0  ,52.0D0  ,52.6D0  ,53.1D0  ,53.6D0  ,54.2D0  ,
     +  54.5D0  ,54.7D0  ,22.6D0   /
      DATA (AE(I, 9, 3),I=1,10)  / 
     +  58.0D0  ,58.4D0  ,59.3D0  ,60.1D0  ,60.7D0  ,61.5D0  ,62.3D0  ,
     +  62.7D0  ,63.1D0  ,22.9D0   /
      DATA (AE(I, 9, 4),I=1,10)  / 
     +  62.0D0  ,63.9D0  ,63.7D0  ,65.7D0  ,65.5D0  ,67.5D0  ,68.2D0  ,
     +  68.9D0  ,69.7D0  ,19.0D0   /
      DATA (AE(I, 9, 5),I=1,10)  / 
     +  .000D+00,72.2D0  ,72.5D0  ,74.2D0  ,74.2D0  ,76.1D0  ,77.0D0  ,
     +  77.8D0  ,78.6D0  ,19.4D0   /
      DATA (AE(I, 9, 6),I=1,10)  / 
     +  .000D+00,80.4D0  ,80.5D0  ,83.1D0  ,83.0D0  ,85.5D0  ,86.8D0  ,
     +  88.1D0  ,89.2D0  ,21.1D0   /
      DATA (AE(I, 9, 7),I=1,10)  / 
     +  .000D+00,63.4D0  ,88.5D0  ,91.3D0  ,91.1D0  ,94.0D0  ,95.8D0  ,
     +  97.3D0  ,98.6D0  ,21.1D0   /
      DATA (AE(I, 9, 8),I=1,10)  / 
     +  .000D+00,.000D+00,98.8D0  ,98.6D0  ,97.8D0  ,102.D0  ,104.D0  ,
     +  106.D0  ,108.D0  ,21.1D0   /
      DATA (AE(I, 9, 9),I=1,10)  / 
     +  .000D+00,.000D+00,84.1D0  ,107.D0  ,107.D0  ,111.D0  ,113.D0  ,
     +  116.D0  ,117.D0  ,22.0D0   /
      DATA (AE(I, 9,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,116.D0  ,115.D0  ,119.D0  ,122.D0  ,
     +  125.D0  ,127.D0  ,23.3D0   /
      DATA (AE(I, 9,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,111.D0  ,123.D0  ,127.D0  ,131.D0  ,
     +  134.D0  ,137.D0  ,24.4D0   /
      DATA (AE(I, 9,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,65.6D0  ,136.D0  ,135.D0  ,140.D0  ,
     +  143.D0  ,146.D0  ,26.3D0   /
      DATA (AE(I, 9,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,146.D0  ,144.D0  ,149.D0  ,
     +  152.D0  ,155.D0  ,27.6D0   /
      DATA (AE(I, 9,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,88.7D0  ,152.D0  ,158.D0  ,
     +  162.D0  ,165.D0  ,29.6D0   /
      DATA (AE(I, 9,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,181.D0  ,167.D0  ,
     +  171.D0  ,174.D0  ,31.8D0   /
      DATA (AE(I, 9,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,117.D0  ,174.D0  ,
     +  180.D0  ,183.D0  ,34.4D0   /
      DATA (AE(I, 9,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,72.0D0  ,201.D0  ,
     +  189.D0  ,192.D0  ,36.4D0   /
      DATA (AE(I, 9,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,151.D0  ,
     +  198.D0  ,201.D0  ,38.3D0   /
      DATA (AE(I, 9,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,95.2D0  ,
     +  220.D0  ,210.D0  ,41.4D0   /
      DATA (AE(I, 9,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  192.D0  ,217.D0  ,43.4D0   /
      DATA (AE(I,10, 1),I=1,10)  / 
     +  62.1D0  ,62.1D0  ,62.6D0  ,62.9D0  ,63.3D0  ,63.3D0  ,64.0D0  ,
     +  64.0D0  ,64.0D0  ,19.0D0   /
      DATA (AE(I,10, 2),I=1,10)  / 
     +  75.1D0  ,75.4D0  ,76.3D0  ,76.8D0  ,77.6D0  ,77.9D0  ,78.8D0  ,
     +  79.0D0  ,79.3D0  ,22.6D0   /
      DATA (AE(I,10, 3),I=1,10)  / 
     +  87.5D0  ,88.3D0  ,89.4D0  ,90.2D0  ,91.3D0  ,91.9D0  ,93.0D0  ,
     +  93.5D0  ,93.9D0  ,22.9D0   /
      DATA (AE(I,10, 4),I=1,10)  / 
     +  104.D0  ,104.D0  ,105.D0  ,106.D0  ,107.D0  ,108.D0  ,109.D0  ,
     +  110.D0  ,110.D0  ,19.0D0   /
      DATA (AE(I,10, 5),I=1,10)  / 
     +  .000D+00,122.D0  ,122.D0  ,123.D0  ,124.D0  ,125.D0  ,126.D0  ,
     +  127.D0  ,128.D0  ,19.4D0   /
      DATA (AE(I,10, 6),I=1,10)  / 
     +  .000D+00,138.D0  ,139.D0  ,140.D0  ,142.D0  ,143.D0  ,144.D0  ,
     +  146.D0  ,147.D0  ,21.1D0   /
      DATA (AE(I,10, 7),I=1,10)  / 
     +  .000D+00,85.3D0  ,158.D0  ,159.D0  ,161.D0  ,162.D0  ,164.D0  ,
     +  166.D0  ,167.D0  ,21.1D0   /
      DATA (AE(I,10, 8),I=1,10)  / 
     +  .000D+00,.000D+00,176.D0  ,177.D0  ,179.D0  ,181.D0  ,183.D0  ,
     +  184.D0  ,186.D0  ,21.1D0   /
      DATA (AE(I,10, 9),I=1,10)  / 
     +  .000D+00,.000D+00,114.D0  ,199.D0  ,201.D0  ,202.D0  ,205.D0  ,
     +  206.D0  ,207.D0  ,22.0D0   /
      DATA (AE(I,10,10),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,218.D0  ,219.D0  ,220.D0  ,224.D0  ,
     +  225.D0  ,226.D0  ,23.3D0   /
      DATA (AE(I,10,11),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,150.D0  ,238.D0  ,238.D0  ,243.D0  ,
     +  244.D0  ,245.D0  ,24.4D0   /
      DATA (AE(I,10,12),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,85.8D0  ,255.D0  ,255.D0  ,261.D0  ,
     +  262.D0  ,263.D0  ,26.3D0   /
      DATA (AE(I,10,13),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,195.D0  ,272.D0  ,279.D0  ,
     +  279.D0  ,280.D0  ,27.6D0   /
      DATA (AE(I,10,14),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,115.D0  ,290.D0  ,296.D0  ,
     +  297.D0  ,298.D0  ,29.6D0   /
      DATA (AE(I,10,15),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,263.D0  ,313.D0  ,
     +  314.D0  ,315.D0  ,31.8D0   /
      DATA (AE(I,10,16),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,150.D0  ,330.D0  ,
     +  331.D0  ,332.D0  ,34.4D0   /
      DATA (AE(I,10,17),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,90.0D0  ,319.D0  ,
     +  349.D0  ,349.D0  ,36.4D0   /
      DATA (AE(I,10,18),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,196.D0  ,
     +  366.D0  ,367.D0  ,38.3D0   /
      DATA (AE(I,10,19),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,122.D0  ,
     +  387.D0  ,384.D0  ,41.4D0   /
      DATA (AE(I,10,20),I=1,10)  / 
     +  .000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,.000D+00,
     +  247.D0  ,401.D0  ,43.4D0     /
      DATA (ERES(I, 1),I=1,10)  / 10*0.D0/
      DATA (ERES(I, 2),I=1,10)  / 10*0.D0/
      DATA (ERES(I, 3),I=1,10)  / 10*0.D0/
      DATA (ERES(I, 4),I=1,10)  / 10*0.D0/
      DATA (ERES(I, 5),I=1,10)  / 10*0.D0/
      DATA (ERES(I, 6),I=1,10)  / 
     +    0.000D0, 0.000D0, 0.000D0, 0.000D0, 0.000D0, 0.000D0, 0.000D0,
     +    2.780D0, 2.880D0, 2.890D0 /
      DATA (ERES(I, 7),I=1,10)  / 
     +    1.500D0, 2.460D0, 2.510D0, 2.610D0, 2.700D0, 2.920D0, 3.070D0,
     +    3.200D0, 3.330D0, 2.890D0 /
      DATA (ERES(I, 8),I=1,10)  / 
     +    4.470D0, 4.350D0, 4.390D0, 4.550D0, 4.660D0, 4.890D0, 4.980D0,
     +    5.100D0, 5.220D0, 2.890D0 /
      DATA (ERES(I, 9),I=1,10)  / 
     +    7.480D0, 7.380D0, 7.370D0, 7.480D0, 7.510D0, 7.630D0, 7.660D0,
     +    7.750D0, 7.820D0, 2.890D0 /
      DATA (ERES(I,10),I=1,10)  / 
     +   15.270D0,15.190D0,15.200D0,15.370D0,15.380D0,15.430D0,15.540D0,
     +   15.590D0,15.630D0, 2.890D0 /
      END
C->
C=======================================================================

      SUBROUTINE FRAGM (IAT,IAP, NW,B, NF, IAF)

C-----------------------------------------------------------------------
C...Nuclear Fragmentation, Abrasion-ablation model, 
C...Based on Jon Engel's routines ABRABL 
C...This most recent version adds for all prefragment
C...masses > 10 the model calculation for the fragment
C...mass distribution and the energy carried by the fragment
C...of W. Friedmann
C...The average values are used to implement the model
C...in the montecarlo fashion / TSS, Dec '91
C.
C.  INPUT: IAP = mass of incident nucleus
C.         IAT = mass of target   nucleus
C.         NW = number of wounded nucleons in the beam nucleus
C.         B  = impact parameter in the interaction
C.     
C.  OUTPUT : NF = number of fragments  of the spectator nucleus
C.           IAF(1:NF) = mass number of each fragment
C.           PF(3,60) in common block /FRAGMENTS/ contains
C.           the three momentum components (MeV/c) of each
C.           fragment in the projectile frame
C..............................................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      COMMON /FRAGMENTS/ PPP(3,60)
      COMMON /FRAGMOD/A(10,10,20),AE(10,10,20),ERES(10,10),NFLAGG(10,10)
      DIMENSION IAF(60)
      DIMENSION AA(10), EAA(10) 
      SAVE
      EXTERNAL GASDEV
      DATA AA/10.D0,15.D0,20.D0,25.D0,30.D0,35.D0,40.D0,45.D0,50.D0,
     $        56.D0/
      DATA EAA/1.D0,2.D0,4.D0,6.D0,8.D0,10.D0,12.D0,16.D0,20.D0,30.D0/

      AP=IAP
      AT=IAT
      NPF = IAP - NW
      IF (NPF .EQ. 0) THEN
         NF = 0
         RETURN
      ENDIF

      EB = ESTAR(AP,AT, B)
      EBP = ESTARP (NPF, NW)
C CONTRIBUTION TO E* FROM ENERGY DEPOSITED BY SECONDARIES
      EB = EB + EBP
C TOTAL E* IS THE SUM OF THE TWO COMPONENTS

C.....Prefragment transverse momentum (MeV/nucleon)...
            FK = FERMK(AP)
C FERMI MOMENTUM OF THE PROJECTILE NUCLEUS
            IF (NW .LT. IAP) THEN
            SIG = FK*DSQRT(NW*NPF/(AP-1.D0))/3.162D0
C GAUSSIAN SIGMA IN ALL THREE DIRECTION
            ELSE
            SIG = FK/3.162D0
C THIS IS NOT CORRECT, TOO LARGE !!!!!!!!!!!!!!
            ENDIF
             PPFX = SIG*GASDEV(0)/NPF
             PPFY = SIG*GASDEV(1)/NPF
C THREE MOMENTUM COMPONENTS PER NUCLEON FOR THE PREFRAGMENT

C.............Crude model for small prefragment mass .......
            IF (NPF .LT. 10) THEN
                 CALL EVAP(NPF, EB, EPS, NNUC, NALP)
C   EPS IS THE KINETIC ENERGY CARRIED BY THE EVAPORATED NUCLEONS
               ETOT = 938.D0 + EPS
                 PP = SQRT((ETOT*ETOT - 8.79844D5)/3.D0)
C   AVERAGE MOMENTUM OF EVAPORATED NUCLEONS IN EACH DIRECTION
                 NUC = NPF - NNUC - 4*NALP
                 NF = 0
                 IF (NUC .GT. 0) THEN
                    NF = NF + 1
                    IAF(NF) = NUC
                    PPP(1,NF) = NUC*PPFX
                    PPP(2,NF) = NUC*PPFY
                 ENDIF
                 IF (NALP .NE. 0) THEN
                 DO I=1,NALP
                   NF = NF + 1
                    IAF(NF) = 4
                   CALL SINCO(S1,C1)
                   CALL SINCO(S2,C2)
                   PXE = 4.D0*PP*S1*S2
                   PYE = 4.D0*PP*S1*C2
                   PPP(1,NF) = 4.D0*PPFX + PXE
                   PPP(2,NF) = 4.D0*PPFY + PYE
                   PPP(1,1) = PPP(1,1) - PXE
                   PPP(2,1) = PPP(2,1) - PYE
                 ENDDO
                 ENDIF
                 IF (NNUC .NE. 0) THEN
                 DO I=1,NNUC
                    NF = NF + 1
                    IAF(NF) = 1
                    CALL SINCO(S1,C1)
                    CALL SINCO(S2,C2)
                    PXE = PP*S1*S2
                    PYE = PP*S1*C2
                    PPP(1,NF) = 4.D0*PPFX + PXE
                    PPP(2,NF) = 4.D0*PPFY + PYE
                    PPP(1,1) = PPP(1,1) - PXE
                    PPP(2,1) = PPP(2,1) - PYE
                 ENDDO
                 ENDIF
                 RETURN
            ENDIF

C.........More refined model calculation .............
      JA = NPF/5 -1
      IF (JA .LT. 10) THEN
      IF ((NPF - AA(JA)) .GT. (AA(JA+1)-NPF)) JA = JA + 1
      ENDIF
      ARAT = DBLE(NPF)/AA(JA)
      DO J=1,10
      IF (EB .LT. EAA(J)) GO TO 29
      ENDDO
      JE = 10
      GO TO 39
   29      JE = J
   39      IF (JE .GT. 1 .AND. JE .NE. 10) THEN
      IF ((EB - EAA(J-1)) .LT. (EAA(J)-EB)) JE = J - 1
      ENDIF
      ERAT = EB/EAA(JE)
        IF (EB .LT. 1.D0) THEN
        ERAT = EB
        ENDIF
C INTERPOLATE BETWEEN EB=0. (NOTHING HAPPENS) AND EB = 1. MeV
         IF (JA .EQ. 10 .AND. JE .GT. 6) THEN
            WRITE(*,*)' JA=',JA,',   JE=',JE
         ENDIF
   43    ESUM = 0.D0
      NSUM = 0
      JF = 0
      DO J=20,1,-1
        FR =  A(JA, JE, J)*ARAT*ERAT
        N1 = INT(1.D0 + FR)
        FR1 = FR/DBLE(N1)
        DO K=1, N1
          IF (S_RNDM(0) .LT. FR1) THEN
            JF = JF + 1
            IAF(JF) = J
            NSUM = NSUM + J
            EKIN = ERAT*AE(JA,JE, J)
            IF (EKIN .GT. 0.D0) THEN
              ESUM = ESUM + EKIN
              ETOT = 938.D0*IAF(JF) + EKIN
              PP = DSQRT(2.D0*(ETOT*ETOT - IAF(JF)**2*8.79844D5)/3.D0)
              CALL SINCO(S1,C1)
              CALL SINCO(S2,C2)
              PPP(1,JF) = PP*S1*S2 + IAF(JF)*PPFX
              PPP(2,JF) = PP*S1*C2 + IAF(JF)*PPFY
            ENDIF
            IF (NSUM .GT. NPF) THEN
C           WRITE(*,*)' WARNING, NSUM=', NSUM,',  NPF=',NPF
C           WRITE(*,*)'  ARAT =', ARAT
              GO TO 43
            ELSE
              IF (NSUM .EQ. NPF) THEN
                GO TO 44
              ENDIF
            ENDIF
          ENDIF
        ENDDO
      ENDDO
      IF (NFLAGG(JA,JE) .EQ. 0) THEN
C 'THE RESIDUE' IS A NUCLEAR FRAGMENT
        JF = JF + 1
        IAF(JF) = NPF - NSUM
        F1 = NPF*EB - ESUM
        IF (F1 .LT. 0.D0) F1 = 0.D0
C GIVE THE REST OF EB TO THE FRAGMENT
        EKIN = F1
        IF (EKIN .GT. 0.D0) THEN
          ETOT = 938.D0*IAF(JF) + EKIN
          PP = DSQRT(2.D0*(ETOT*ETOT - IAF(JF)**2*8.79844D5)/3.D0)
          CALL SINCO(S1,C1)
          CALL SINCO(S2,C2)
          PPP(1,JF) = PP*S1*S2 + IAF(JF)*PPFX
          PPP(2,JF) = PP*S1*C2 + IAF(JF)*PPFY
        ENDIF
      ELSE
C 'THE RESIDUE' CONSISTS OF SPECTATOR NUCLEONS
        N1 = NPF - NSUM
        DO K=1,N1
          JF = JF + 1
          IAF(JF) = 1
          EKIN = ERAT*ERES(JA,JE)
          IF (EKIN .GT. 0.D0) THEN
            ETOT = 938.D0*IAF(JF) + EKIN
            PP = DSQRT(2.D0*(ETOT*ETOT - IAF(JF)**2*8.79844D5)/3.D0)
            CALL SINCO(S1,C1)
            CALL SINCO(S2,C2)
            PPP(1,JF) = PP*S1*S2 + PPFX
            PPP(2,JF) = PP*S1*C2 + PPFY
          ENDIF
        ENDDO
      ENDIF
  44  NF = JF
      RETURN
      END
C->
C=======================================================================

      FUNCTION ESTARP (NPF, NW)

C-----------------------------------------------------------------------
C CONTRIBUTION TO E* FROM ENERGY DEPOSITED BY SECONDARIES
C VERY NAIVE VERSION INCORPORATING HUEFFNER'S IDEAS
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      SAVE

      APF = NPF
      F1 = 15.3D0/APF**0.666666666D0
C AVERAGE KINETIC ENERGY/NUCLEON IN PREFRAGMENT (MeV)
C PER PATHLENGTH EQUAL TO THE PREFRAGMENT RADIUS
      ESTARP = 0.D0
      DO I=1,NW
        IF (S_RNDM(0) .GT. 0.5D0) THEN
          F2 = F1*RDIS(0)
          ESTARP = ESTARP + F2
        ENDIF
      ENDDO
C SAMPLE RANDOMLY PER WOUNDED NUCLEON, x NW
      RETURN
      END
C=======================================================================
      
      FUNCTION RDIS(Idum)

C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      dimension probr(20)
      SAVE
      data probr/
     *      0.10000D0, 0.15748D0, 0.21778D0, 0.28605D0, 0.36060D0,
     *      0.43815D0, 0.51892D0, 0.60631D0, 0.70002D0, 0.79325D0,
     *      0.88863D0, 0.98686D0, 1.10129D0, 1.21202D0, 1.32932D0,
     *      1.44890D0, 1.57048D0, 1.70139D0, 1.83417D0, 2.00000D0/

      rdis = idum
      nr = INT(20.D0*S_RNDM(0) + 1.D0)
      if (nr .eq. 1) then
        f1 = 0.D0
      else
        f1 = probr(nr-1)
      endif
      dr = probr(nr) - f1
      rdis = f1 + dr*S_RNDM(1)
      return
      end

C=======================================================================

      FUNCTION ESTAR(ap,at,b)

C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN

      SAVE

c      real*4 ap,at,b,estar
      sigma=4.5D0  !total n-n cross section in fm**2
      rt=.82d0*at**.33333333D0 !target radius
      rp=.82d0*ap**.33333333D0 !projectile radius
      alpha=rt**2/rp**2
      beta=b**2/rt**2
      f=at*sigma/(PI*rt**2)
      alf = log(f)
      alalf = log(alpha)
      gfac=0.d0
      gfac1=0.d0
      s1=0.D0
      s2=0.D0
      s3=0.D0
      ii=1
      do n=0,10 ! This limit may not need to be so high.
         if(n.ge.2) then
            gfac1=gfac
            gfac=gfac+log(float(n)) 
         endif
         g0=n*alf -n*beta*alpha/(n+alpha)+alalf
         g1=g0-log(alpha+n)-gfac
         g2=(n+2)*log(f)-(n+2)*beta*alpha/(n+2+alpha) 
     >      +log(n+2+alpha+beta*alpha**2)-3.d0*log(n+2.d0+alpha)-gfac
         g3=g0-2.d0*log(n+alpha)-gfac1
         ii=-ii
         s1=s1+ii*exp(g1)
         s2=s2+ii*exp(g2)
         if(n.ge.1) s3=s3+ii*exp(g3)
      enddo

      pb=s1
      e1b=197.D0**2/(2.D0*938.d0*rp**2*pb) *s2
c      a=b*(s3/pb-1)
c      a=-b*s3/pb
c      e2b=-.5* 938. * (41./(ap**.333))**2 * a**2 /(197.**2)
c      estar=e1b+e2b
      estar = e1b
      return
      end
C=======================================================================

      SUBROUTINE EVAP(npf,eb,eps,nnuc,nalp)

C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      SAVE

      eps=7.5D0+sqrt(8.D0*eb)
      n=min(npf*int(eb/eps),npf)
      nalp=n/5
      nnuc=n-4*nalp
      return
      end
C->
C=======================================================================

      FUNCTION FERMK(A)

C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      DIMENSION AA(6), FK(6)
      SAVE
      DATA AA/4.D0, 6.D0, 12.D0, 24.D0, 40.D0, 57.D0/
      DATA FK/130.D0,169.D0,221.D0,235.D0,251.D0,260.D0/

      DO I=2,4
      IF (A .LT. AA(I)) GO TO 25
      ENDDO
      I = 5
   25      F11 = AA(I-1)
      F12 = AA(I)
      F13 = AA(I+1)
      F21 = FK(I-1)
      F22 = FK(I)
      F23 = FK(I+1)
      FERMK = QUAD_INT(A,F11,F12,F13, F21,F22,F23)
      RETURN
      END

C*=======================================================================
C. Multiple interaction structure
C========================================================================

      SUBROUTINE INT_NUC (IA, IB, SIG0, SIGEL) 

C-----------------------------------------------------------------------
C...Compute with a montecarlo code  the  "multiple interaction structure"
C.  of a nucleus-nucleus interaction
C.
C.  INPUT : IA            = mass of target nucleus
C.          IB            = mass of projectile nucleus
C.          SIG0 (mbarn)  = inelastic pp cross section
C.          SIGEL(mbarn)  = elastic pp cross section
C.
C.  OUTPUT : in common block /CNUCMS/
C.           B = impact parameter (fm)
C.           BMAX = maximum impact parameter for generation
C.           NTRY = number of "trials" before one interaction
C.           NA = number of wounded nucleons in A
C.           NB =    "        "        "     in B
C.           NI = number of nucleon-nucleon inelastic interactions 
C.           NAEL = number of elastically scattered nucleons in  A 
C.           NBEL =    "         "           "          "    in  B
C.           JJA(J)  [J=1:IA]   = number of inelastic interactions 
C.                                of J-th nucleon of nucleus A
C.           JJB(J)  [J=1:IB]   = number of inelastic interactions 
C.                                of J-th nucleon of nucleus B
C.           JJAEL(J)  [J=1:IA]   = number of elastic interactions 
C.                                of J-th nucleon of nucleus A
C.           JJBEL(J)  [J=1:IB]   = number of elastic interactions 
C.                                of J-th nucleon of nucleus B
C.           JJINT(J,K)  [J=1:NB, K=1:NA]  (0 = no interaction) 
C.                                         (1 = interaction )
C.                                         between nucleon J of A and K of B
C-----------------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      PARAMETER (IAMAX=56)
      COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
     +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
     +         ,JJAEL(IAMAX), JJBEL(IAMAX)
      DIMENSION XA(IAMAX), YA(IAMAX), XB(IAMAX), YB(IAMAX)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

      SIGT = SIG0 + SIGEL
      R2  = 0.1D0 * SIG0/PI
      R2T = 0.1D0 * SIGT/PI
      BMAX = 15.D0                             ! fm
      NTRY = 0
      CALL NUC_CONF (IA, XA, YA)
      CALL NUC_CONF (IB, XB, YB)
      NI = 0
      NIEL = 0
      DO JA=1,IA
         JJA(JA) = 0
         JJAEL(JA) = 0
      ENDDO
      DO JB=1,IB
         JJB(JB) = 0
         JJBEL(JB) = 0
         DO JA=1,IA
            JJINT(JB,JA) = 0
         ENDDO
      ENDDO
1000  B = BMAX*SQRT(S_RNDM(0))
      PHI = TWOPI*S_RNDM(1)
      BX = B*COS(PHI)
      BY = B*SIN(PHI)
      NTRY = NTRY+1
      DO JA=1,IA
         DO JB=1,IB
            S = (XA(JA)-XB(JB)-BX)**2 + (YA(JA)-YB(JB)-BY)**2
            IF (S .LT. R2)  THEN
               NI = NI + 1
               JJA(JA) = JJA(JA)+1
               JJB(JB) = JJB(JB)+1
               JJINT(JB,JA) = 1
            ELSE IF (S .LT. R2T)  THEN
               NIEL = NIEL + 1
               JJAEL(JA) = JJAEL(JA)+1
               JJBEL(JB) = JJBEL(JB)+1
            ENDIF
         ENDDO
      ENDDO
      IF (NI + NIEL .EQ. 0)  GOTO 1000
      NA = 0
      NB = 0
      NAEL = 0
      NBEL = 0
      DO JA=1,IA
         IF (JJA(JA) .GT. 0)  THEN
            NA = NA + 1
         ELSE
            IF (JJAEL(JA) .GT. 0)  NAEL = NAEL+1
         ENDIF
      ENDDO
      DO JB=1,IB
         IF (JJB(JB) .GT. 0)  THEN
            NB = NB + 1
         ELSE
            IF (JJBEL(JB) .GT. 0)  NBEL = NBEL+1
         ENDIF
      ENDDO
      RETURN
      END
C=======================================================================

       SUBROUTINE NUC_CONF (IA, XX, YY)

C-----------------------------------------------------------------------
C...This routine generates the configuration  of a nucleus 
C.  need an initialization call to NUC_GEOM_INI
C.
C.  INPUT  : IA = mass number of the nucleus
C.  OUTPUT : XX(1:IA), YY(1:IA) (fm) = position in impact parameter
C.                                     space of the IA nucleons
C...................................................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      PARAMETER (IAMAX=56)
      DIMENSION XX(IAMAX), YY(IAMAX)
      PARAMETER (NB=401)
      COMMON /CPROFA/ ZMIN, DZ, BBZ(NB,IAMAX)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

      DO J=1,IA
         Z = S_RNDM(J)
         JZ = INT((Z-ZMIN)/DZ)+1
         JZ = MIN(JZ,400)
         T = (Z-ZMIN)/DZ - DBLE(JZ-1)
         B = BBZ(JZ,IA)*(1.D0-T) + BBZ(JZ+1,IA)*T
         PHI = TWOPI*S_RNDM(J+1)
         XX(J) = B*COS(PHI)
         YY(J) = B*SIN(PHI)
      ENDDO
      RETURN
      END
C=======================================================================

      SUBROUTINE NUC_GEOM_INI

C-----------------------------------------------------------------------
C...Initialize all nucleus profiles
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      PARAMETER (NB=401)
      PARAMETER (IAMAX=56)
      COMMON /CPROF/ DB, BMAX, BB(NB), TB(NB), A
      COMMON /CPROFA/ ZMIN, DZ, BBZ(NB,IAMAX)
      DIMENSION FFB(NB), GGB(NB)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

      CALL SHELL_INI
      CALL WOOD_SAXON_INI
      DO IA= 2,IAMAX
           JA = IA
         CALL NUC_PROFIL(JA)
         DO K=1,NB
           FFB(K) = BB(K)*TB(K) * TWOPI
         ENDDO            
         GGB(1) = 0.D0
         GGB(NB) = 1.D0
         DO K=2,NB-1
           GGB(K) = GGB(K-1) + FFB(K-1)*DB
         ENDDO            
         CALL INVERT_ARRAY(GGB,0.D0,DB,NB, BBZ(1,IA), ZMIN, DZ)
      ENDDO
      RETURN
      END
C=======================================================================

      SUBROUTINE NUC_PROFIL (JA)

C-----------------------------------------------------------------------
C...Compute the profile function T(b)
C.  normalised as INT[d2b T(b) = 1]
C.  INPUT : JA = integer mass number of nucleus
C...............................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      PARAMETER (NB=401)
      EXTERNAL DENSA
      DOUBLE PRECISION DENSA
      COMMON /CC01/  B
      COMMON /CCDA/ JJA
      COMMON /CPROF/ DB, BMAX, BB(NB), TB(NB), A
      SAVE

      BMAX = 7.5D0
      DB = BMAX/DBLE(NB-1)
      JJA = JA
      A = JA
      DO JB=1,NB
        B = DB*DBLE(JB-1)
        BB(JB) = B
        IF (JA .LE. 18)  THEN
            TB(JB) = PROFNUC (B, JA)
         ELSE
            TB(JB) = 2.D0*GAUSS (DENSA,0.D0,BMAX)
         ENDIF
      ENDDO
      RETURN
      END
C=======================================================================

      SUBROUTINE NUC1_PROFIL (AA)

C-----------------------------------------------------------------------
C...Compute the profile function T(b)
C.  normalised as INT[d2b T(b) = 1]
C.  INPUT : AA = mass number of nucleus
C...............................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      PARAMETER (NB=401)
      EXTERNAL DENSA
      DOUBLE PRECISION DENSA
      COMMON /CC01/  B
      COMMON /CPROF/ DB, BMAX, BB(NB), TB(NB), A
      SAVE

      A = AA
      IA1 = INT(AA)
      IA2 = IA1 + 1
      U = AA - DBLE(IA1)
      BMAX = 7.5D0
      DB = BMAX/DBLE(NB-1)
      DO JB=1,NB
         B = DB*DBLE(JB-1)
         BB(JB) = B
         IF (A .LE. 18.D0)  THEN
             T1 = PROFNUC (B, IA1)
             T2 = PROFNUC (B, IA2)
          ELSE
             JJA = IA1
             T1 = 2.D0*GAUSS (DENSA,0.D0,BMAX)
             JJA = IA2
             T2 = 2.D0*GAUSS (DENSA,0.D0,BMAX)
          ENDIF
          TB(JB) = (1.D0-U)*T1  + U*T2
      ENDDO
      RETURN
      END

C*======================================================================
C.   Code about nuclear densities
C=======================================================================

      FUNCTION DENS_NUC (R, JA)

C-----------------------------------------------------------------------
C....Nuclear density (normalised to 1)
C.   for a nucleus of mass number JA
C.   INPUT R = radial coordinate  (fm)
C.         JA = integer mass number
C.  OUTPUT (fm**-3)
C--------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
      SAVE

      IF (JA .GT. 18)  THEN
         DENS_NUC = WOOD_SAXON(R,JA)
      ELSE IF (JA .NE. 4)  THEN
         DENS_NUC = HELIUM(R)
      ELSE
         DENS_NUC = SHELL(R,JA)
      ENDIF
      RETURN
      END
C=======================================================================

      FUNCTION WOOD_SAXON (R, JA) 

C-----------------------------------------------------------------------
C....Wood-Saxon nuclear density (normalised to 1)
C.   for a nucleus of mass number A.
C.   INPUT R =  (fm)
C.         JA = mass number
C.   OUTPUT (fm**-3)
C------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
      SAVE

      WOOD_SAXON = CC0(JA)/(1.D0+EXP((R-RR0(JA))/AA0(JA)))
      RETURN
      END      
C=======================================================================

      FUNCTION HELIUM (R)

C-----------------------------------------------------------------------
C... Helium density from Barrett and Jackson
C.   INPUT R = r coordinate (fm)
C.   OUTPUT (fm**-3)
C........................................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      SAVE
      DATA R0 /0.964D0/, CA /0.322D0/   ! fm
      DATA W /0.517D0/, CC /5.993224D-02/

      HELIUM = CC*(1.D0+W*(R/R0)**2)/(1.D0 + EXP((R-R0)/CA))
      RETURN
      END
C=======================================================================

      FUNCTION SHELL (R,JA)

C-----------------------------------------------------------------------
C...Density in the shell model
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON /CSHELL/ RR0(18), RR02(18)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

      R0 = RR0(JA)
      C1 = MIN(1.D0,4.D0/DBLE(JA))
      CS = 1.D0/(R0**3 * PI**1.5D0)
      CP = 2.D0*CS/3.D0
      FS = EXP(-(R/R0)**2)
      FP = (R/R0)**2 * FS
      SHELL = C1*CS*FS + (1.D0-C1)*CP*FP
      RETURN
      END
C=======================================================================

      FUNCTION PROFNUC (B, JA)

C-----------------------------------------------------------------------
C...This function return
C.  the profile T(b) for a nucleus of mass number A
C.  INPUT B = impact parameter (GeV**-1)
C.        JA = integer mass number
C.  OUTPUT  (fm**-2)
C.
C.  The  density of the nucleus is the `shell model density'
C.  the parameter r0 must beinitialized in the common block
C.............................................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON /CSHELL/ RR0(18), RR02(18)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

      B2 = B*B
      ARG = B2/RR02(JA)
      TS = EXP(-ARG)
      TP = TS*(2.D0*B2+RR02(JA))/(3.D0*RR02(JA))
      CS = MIN(1.D0,4.D0/DBLE(JA))
      PROFNUC = (CS*TS + (1.D0-CS)*TP)/(PI*RR02(JA))
      RETURN
      END
C=======================================================================

      SUBROUTINE SHELL_INI

C-----------------------------------------------------------------------
C...Initialize the parameter  of the shell model
C.  for the nuclei with    6 < A < 18
C..............................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      COMMON /CSHELL/ RR0(18), RR02(18)
      DIMENSION RR(18)
      SAVE
C...Data on Sqrt[<r**2>]  in fermi
      DATA RR /0.81D0, 2.095D0,  1.88D0, 1.674D0,  -1.D0,
     +         2.56D0, 2.41D0,    -1.D0, 2.519D0, 2.45D0,
     +         2.37D0, 2.460D0, 2.440D0,  2.54D0, 2.58D0, 
     +         2.718D0,2.662D0, 2.789D0/

      DO JA=1,18
         A = DBLE(JA)
         RMED = RR(JA)
         IF (RMED .LE. 0.D0)   RMED = 0.5D0*(RR(JA-1) + RR(JA+1))
         C = MAX(1.5D0,(5.D0/2.D0 - 4.D0/A) )
         R0 = RMED/SQRT(C)
         RR0 (JA) = R0
         RR02(JA) = R0*R0
      ENDDO
      RETURN
      END
C->
C=======================================================================

      SUBROUTINE WOOD_SAXON_INI

C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

C...Wood-Saxon parameters from  table 6.2   of Barrett and Jackson
      RR0 (19) = 2.59D0
      AA0 (19) = 0.564D0
      RR0 (20) = 2.74D0
      AA0 (20) = 0.569D0
      RR0 (22) = 2.782D0
      AA0 (22) = 0.549D0
      RR0 (24) = 2.99D0
      AA0 (24) = 0.548D0
      RR0 (27) = 2.84D0
      AA0 (27) = 0.569D0
      RR0 (28) = 3.14D0
      AA0 (28) = 0.537D0
      RR0 (29) = 3.77D0
      AA0 (29) = 0.52D0
      RR0 (48) = 3.912D0
      AA0 (48) = 0.5234D0
      RR0 (56) = 3.98D0
      AA0 (56) = 0.569D0
      DO J=19, 56
         IF (RR0(J) .LE. 0.D0)  THEN
            RR0(J) = 1.05D0*DBLE(J)**0.333333333333D0
            AA0(J) = 0.545D0
         ENDIF
         CC0(J)=3.D0/(4.D0*PI*RR0(J)**3)/(1.D0+((AA0(J)*PI)/RR0(J))**2)
      ENDDO
      RETURN
      END
C=======================================================================

      FUNCTION DENSA (Z)

C-----------------------------------------------------------------------
C....Woods Saxon nuclear density (normalised to 1)
C.   for a nucleus of mass number A.
C.   INPUT z = z coordinate (fm)
C.         JA = integer mass number
C.         B (in common /CC01/)  impact parameter  (fm)
C.  OUTPUT (fm**-3)
C--------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON /CC01/  B
      COMMON /CCDA/ JA
      COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
      SAVE

      R = SQRT (Z*Z + B*B)
      DENSA = CC0(JA)/(1.D0+EXP((R-RR0(JA))/AA0(JA)))
      RETURN
      END

C*=====================================================================
C. Cross sections
C======================================================================

      SUBROUTINE SIGMA_AIR (IB,SIG0,SIGEL,KINT,
     +                            SIGMA,DSIGMA,SIGQE,DSIGQE)

C-----------------------------------------------------------------------
C...Compute with a montecarlo method the "production"
C.  and "quasi-elastic" cross section for  
C.  a nucleus-air  interaction 
C.
C.  INPUT : IB            = mass of projectile nucleus
C.          SIG0 (mbarn)  = inelastic pp cross section
C.          KINT            = number  of interactions to generate
C.  OUTPUT : SIGMA (mbarn) = "production" cross section
C.           DSIGMA   "    = error
C.           SIGQE    "    = "quasi-elastic" cross section
C.           DSIGQE   "    = error
C.           additional output is in the common block  /CPROBAB/
C..........................................................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      PARAMETER (IAMAX=56)
      PARAMETER (IAMAX2=3136)          ! IAMAX*IAMAX
      COMMON  /CPROBAB/ PROBA(IAMAX), DPROBA(IAMAX), 
     +   PROBB(IAMAX), DPROBB(IAMAX), PROBI(IAMAX2), DPROBI(IAMAX2),
     +   P1AEL(0:IAMAX),DP1AEL(0:IAMAX),P1BEL(0:IAMAX), DP1BEL(0:IAMAX),
     +   P2AEL(0:IAMAX),DP2AEL(0:IAMAX),P2BEL(0:IAMAX), DP2BEL(0:IAMAX)
      COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
     +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
     +         ,JJAEL(IAMAX), JJBEL(IAMAX)
      DIMENSION  MMA(0:IAMAX), MMB(0:IAMAX), MMI(0:IAMAX2)
      DIMENSION  M1AEL(0:IAMAX), M1BEL(0:IAMAX)
      DIMENSION  M2AEL(0:IAMAX), M2BEL(0:IAMAX)
      DOUBLE PRECISION FOX
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE
      DATA FOX /0.21522D0/  !atomic percentage of 'non-nitrogen' in air

      R2 = 0.1D0 * SIG0/PI
      BMAX = 15.D0                           ! fm
      SIGMA0 = PI*BMAX*BMAX*10.              ! mbarn
      IA = 16
      DO J=1,IA
         MMA(J) = 0
         M1AEL(J) = 0
         M2AEL(J) = 0
      ENDDO
      DO J=1,IB
         MMB(J) = 0
         M1BEL(J) = 0
         M2BEL(J) = 0
      ENDDO
      DO J=1,IA*IB
         MMI(J) = 0
      ENDDO
      NN = 0
      M = 0
      DO KK=1,KINT
c  select target IA from air composition
         R = S_RNDM(KK)
         IA = 14
         IF (R .LT. FOX)  IA = 16

         CALL INT_NUC (IA, IB, SIG0, SIGEL) 
         NN = NN + NTRY
         MMI(NI) = MMI(NI) + 1
         MMA(NA) = MMA(NA)+1
         MMB(NB) = MMB(NB)+1
         IF (NI .GT. 0)  THEN
            M = M+1
            M1AEL(NAEL) = M1AEL(NAEL)+1
            M1BEL(NBEL) = M1BEL(NBEL)+1
         ELSE
            M2AEL(NAEL) = M2AEL(NAEL)+1
            M2BEL(NBEL) = M2BEL(NBEL)+1
         ENDIF
      ENDDO
      MQE = KINT - M
      SIGMA  = SIGMA0 * DBLE(M)/DBLE(NN)
      DSIGMA = SIGMA0 * SQRT(DBLE(M))/DBLE(NN)
      SIGQE  = SIGMA0 * DBLE(MQE)/DBLE(NN)
      DSIGQE = SIGMA0 * SQRT(DBLE(MQE))/DBLE(NN)
      DO J=1,IA
         PROBA(J) = DBLE(MMA(J))/DBLE(M)
         DPROBA(J) = SQRT(DBLE(MMA(J)))/DBLE(M)
      ENDDO
      DO J=1,IB
         PROBB(J) = DBLE(MMB(J))/DBLE(M)
         DPROBB(J) = SQRT(DBLE(MMB(J)))/DBLE(M)
      ENDDO
      DO J=1,IA*IB
         PROBI(J) = DBLE(MMI(J))/DBLE(M)
         DPROBI(J) = SQRT(DBLE(MMI(J)))/DBLE(M)
      ENDDO
      DO J=0,IA
         P1AEL(J) = DBLE(M1AEL(J))/DBLE(M)
         DP1AEL(J) = SQRT(DBLE(M1AEL(J)))/DBLE(M)
         P2AEL(J) = DBLE(M2AEL(J))/DBLE(MQE)
         DP2AEL(J) = SQRT(DBLE(M2AEL(J)))/DBLE(MQE)
      ENDDO
      DO J=0,IB
         P1BEL(J) = DBLE(M1BEL(J))/DBLE(M)
         DP1BEL(J) = SQRT(DBLE(M1BEL(J)))/DBLE(M)
         P2BEL(J) = DBLE(M2BEL(J))/DBLE(MQE)
         DP2BEL(J) = SQRT(DBLE(M2BEL(J)))/DBLE(MQE)
      ENDDO
      RETURN
      END
C->
C=======================================================================

      SUBROUTINE SIGMA_MC (IA,IB,SIG0,SIGEL,KINT,
     +                            SIGMA,DSIGMA,SIGQE,DSIGQE)

C-----------------------------------------------------------------------
C...Compute with a montecarlo method the "production"
C.  and "quasi-elastic" cross section for  
C.  a nucleus-nucleus interaction
C.
C.  INPUT : IA            = mass of target nucleus
C.          IB            = mass of projectile nucleus
C.          SIG0 (mbarn)  = inelastic pp cross section
C.          KINT            = number  of interactions to generate
C.  OUTPUT : SIGMA (mbarn) = "production" cross section
C.           DSIGMA   "    = error
C.           SIGQE    "    = "quasi-elastic" cross section
C.           DSIGQE   "    = error
C.           additional output is in the common block  /CPROBAB/
C.           Prob(n_A), Prob(n_B), Prob(n_int)
C..........................................................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      PARAMETER (IAMAX=56)
      PARAMETER (IAMAX2=3136)          ! IAMAX*IAMAX
      COMMON  /CPROBAB/ PROBA(IAMAX), DPROBA(IAMAX), 
     +   PROBB(IAMAX), DPROBB(IAMAX), PROBI(IAMAX2), DPROBI(IAMAX2),
     +   P1AEL(0:IAMAX),DP1AEL(0:IAMAX),P1BEL(0:IAMAX), DP1BEL(0:IAMAX),
     +   P2AEL(0:IAMAX),DP2AEL(0:IAMAX),P2BEL(0:IAMAX), DP2BEL(0:IAMAX)
      COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
     +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
     +         ,JJAEL(IAMAX), JJBEL(IAMAX)
      DIMENSION  MMA(0:IAMAX), MMB(0:IAMAX), MMI(0:IAMAX2)
      DIMENSION  M1AEL(0:IAMAX), M1BEL(0:IAMAX)
      DIMENSION  M2AEL(0:IAMAX), M2BEL(0:IAMAX)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

      R2 = 0.1D0 * SIG0/PI
      BMAX = 15.D0                           ! fm
      SIGMA0 = PI*BMAX*BMAX*10.D0              ! mbarn
      DO J=1,IA
         MMA(J) = 0
         M1AEL(J) = 0
         M2AEL(J) = 0
      ENDDO
      DO J=1,IB
         MMB(J) = 0
         M1BEL(J) = 0
         M2BEL(J) = 0
      ENDDO
      DO J=1,IA*IB
         MMI(J) = 0
      ENDDO
      NN = 0
      M = 0
      DO KK=1,KINT
         CALL INT_NUC (IA, IB, SIG0, SIGEL) 
         NN = NN + NTRY
         MMI(NI) = MMI(NI) + 1
         MMA(NA) = MMA(NA)+1
         MMB(NB) = MMB(NB)+1
         IF (NI .GT. 0)  THEN
            M = M+1
            M1AEL(NAEL) = M1AEL(NAEL)+1
            M1BEL(NBEL) = M1BEL(NBEL)+1
         ELSE
            M2AEL(NAEL) = M2AEL(NAEL)+1
            M2BEL(NBEL) = M2BEL(NBEL)+1
         ENDIF
      ENDDO
      MQE = KINT - M
      SIGMA  = SIGMA0 * DBLE(M)/DBLE(NN)
      DSIGMA = SIGMA0 * SQRT(DBLE(M))/DBLE(NN)
      SIGQE  = SIGMA0 * DBLE(MQE)/DBLE(NN)
      DSIGQE = SIGMA0 * SQRT(DBLE(MQE))/DBLE(NN)
      DO J=1,IA
         PROBA(J) = DBLE(MMA(J))/DBLE(M)
         DPROBA(J) = SQRT(DBLE(MMA(J)))/DBLE(M)
      ENDDO
      DO J=1,IB
         PROBB(J) = DBLE(MMB(J))/DBLE(M)
         DPROBB(J) = SQRT(DBLE(MMB(J)))/DBLE(M)
      ENDDO
      DO J=1,IA*IB
         PROBI(J) = DBLE(MMI(J))/DBLE(M)
         DPROBI(J) = SQRT(DBLE(MMI(J)))/DBLE(M)
      ENDDO
      DO J=0,IA
         P1AEL(J) = DBLE(M1AEL(J))/DBLE(M)
         DP1AEL(J) = SQRT(DBLE(M1AEL(J)))/DBLE(M)
         P2AEL(J) = DBLE(M2AEL(J))/DBLE(MQE)
         DP2AEL(J) = SQRT(DBLE(M2AEL(J)))/DBLE(MQE)
      ENDDO
      DO J=0,IB
         P1BEL(J) = DBLE(M1BEL(J))/DBLE(M)
         DP1BEL(J) = SQRT(DBLE(M1BEL(J)))/DBLE(M)
         P2BEL(J) = DBLE(M2BEL(J))/DBLE(MQE)
         DP2BEL(J) = SQRT(DBLE(M2BEL(J)))/DBLE(MQE)
      ENDDO
      RETURN
      END

C*=============================================================
C.  Cross sections
C*=============================================================

C Glauber h-air cross section calculation moved to inelScreen src file..

C-----------------------------------------------------------------------
C.  Fit of Block and Cahn to pp and pbar-p cross sections
C-----------------------------------------------------------------------
C=======================================================================

      SUBROUTINE BLOCK(SQS,SIG1,SIG2,SLOP1,SLOP2,
     +                 RHO1,RHO2,SIGEL1,SIGEL2)

C-----------------------------------------------------------------------
C...p-p and pbar-p cross sections
C.  Parametrization of  Block and Cahn
C
C.  INPUT  : SQS   (GeV)  = c.m. energy
C.  
C.  OUPUT : SIG1 (mbarn)    = pp  total  cross section 
C.          SLOP1 (GeV**2)  = slope of elastic scattering
C.          RHO1            = Real/Imaginary part of the amplitude
C.                            for forward elastic  scattering (pp)
C.          SIGEL1 (mbarn)  = pp  elastic scattering  cross section
C.          [1 -> 2   : pp -> pbar p]
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

      S = SQS*SQS
      CALL FPLUS  (S, FR, FI)
      CALL FMINUS (S, GR, GI)
      SIG1 = FI-GI
      SIG2 = FI+GI
      RHO1 = (FR-GR)/(FI-GI)
      RHO2 = (FR+GR)/(FI+GI)
      CALL SSLOPE (S, BP, BM)
      SLOP1 = BP - GI/FI*(BM-BP)
      SLOP2 = BP + GI/FI*(BM-BP)
      SIGEL1 = SIG1**2*(1.D0+RHO1**2)/(16.D0*PI*SLOP1)/CMBARN
      SIGEL2 = SIG2**2*(1.D0+RHO2**2)/(16.D0*PI*SLOP2)/CMBARN
      RETURN
      END
C=======================================================================

      SUBROUTINE FPLUS (S, FR, FI)

C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON /BLOCKC/ AA, BETA, S0, CC, AMU, DD, ALPHA, A0
      COMPLEX*16 Z1, Z2, Z3
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

      F1 = LOG(S/S0)
      Z1 = DCMPLX(F1,-PI/2.D0)
      Z1 = Z1*Z1
      Z2 = 1.D0 + A0*Z1
      Z3 = Z1/Z2
      F2 = CC*S**(AMU-1.D0)
      F3 = 0.5D0*PI*(1.-AMU)
      FI = AA + F2*COS(F3) + BETA*DREAL(Z3)
      FR = -BETA*DIMAG(Z3)+F2*SIN(F3)
      RETURN
      END
C=======================================================================

      SUBROUTINE FMINUS (S, FR, FI)

C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON /BLOCKC/ AA, BETA, S0, CC, AMU, DD, ALPHA, A0
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

      F1 = S**(ALPHA-1.D0)
      F2 = 0.5D0*PI*(1.D0-ALPHA)
      FR = -DD*F1*COS(F2)
      FI = -DD*F1*SIN(F2)
      RETURN
      END
C=======================================================================

      SUBROUTINE SSLOPE (S, BP, BM)

C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON /BLOCKD/ CP, DP, EP, CM, DM
      SAVE

      AL = LOG(S)
      BP = CP + DP*AL + EP*AL*AL
      BM = CM + DM*AL
      RETURN
      END
C=======================================================================

      SUBROUTINE BLOCK_INI

C-----------------------------------------------------------------------
C...Parameters of fit IFIT=1 of Block and Cahn
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON /BLOCKC/ AA, BETA, S0, CC, AMU, DD, ALPHA, A0
      COMMON /BLOCKD/ CP, DP, EP, CM, DM
      SAVE

      AA = 41.74D0
      BETA = 0.66D0
      S0 = 338.5D0
      CC = 0.D0
      AMU = 0.D0
      DD = -39.37D0
      ALPHA = 0.48D0
      A0 = 0.D0
      CP = 10.90D0
      DP = -0.08D0
      EP = 0.043D0
      CM = 23.27D0
      DM = 0.93D0
      RETURN
      END

C*=============================================================
C.  Nucleus-nucleus cross sections
C=======================================================================

      SUBROUTINE SIGNUC_INI (IA,E0)

C-----------------------------------------------------------------------
C...This subroutine receives in INPUT E0 (TeV)
C.  energy per nucleon and computes the cross sections
C.  and interactions lengths for  all nuclei
C.  with A  between 2 and IA
C.  The output is contained in common block /CLENNN/
C.
C.  Attention: the tabulated cross sections are obtained with
C.  new p-p cross sections as used in SIBYLL 2x,
C.  in addition field dimensions changed (RE 04/2000)
C.
C........................................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON /CLENNN/ SSIGNUC(60), ALNUC(60)
      DIMENSION SIGMA(6,56), SIGQE(6,56)
      DIMENSION AA(6)
      SAVE
      DATA NE /6/, AMIN /1.D0/, DA /1.D0/
      DATA AA /1.D0,2.D0,3.D0,4.D0,5.D0,6.D0/
      DATA AVOG /6.0221367D-04/
      DATA ATARGET /14.514D0/            ! effective masss of air
C...Data on `inelastic-production' nucleus-air cross section
      DATA (SIGMA(J, 2),J=1,6) /
     &3.842D+02,4.287D+02,4.940D+02,5.887D+02,6.922D+02,7.767D+02/
      DATA (SIGMA(J, 3),J=1,6) /
     &4.601D+02,5.149D+02,5.595D+02,6.663D+02,7.641D+02,8.446D+02/
      DATA (SIGMA(J, 4),J=1,6) /
     &4.881D+02,5.373D+02,6.005D+02,6.895D+02,7.716D+02,8.967D+02/
      DATA (SIGMA(J, 5),J=1,6) /
     &5.874D+02,6.176D+02,7.181D+02,7.993D+02,9.089D+02,1.031D+03/
      DATA (SIGMA(J, 6),J=1,6) /
     &7.054D+02,7.399D+02,8.388D+02,9.463D+02,1.080D+03,1.197D+03/
      DATA (SIGMA(J, 7),J=1,6) /
     &7.192D+02,7.611D+02,8.449D+02,9.539D+02,1.061D+03,1.176D+03/
      DATA (SIGMA(J, 8),J=1,6) /
     &7.550D+02,7.975D+02,9.153D+02,9.944D+02,1.126D+03,1.236D+03/
      DATA (SIGMA(J, 9),J=1,6) /
     &7.929D+02,8.392D+02,9.265D+02,1.059D+03,1.167D+03,1.262D+03/
      DATA (SIGMA(J, 10),J=1,6) /
     &8.157D+02,8.644D+02,9.512D+02,1.058D+03,1.182D+03,1.298D+03/
      DATA (SIGMA(J, 11),J=1,6) /
     &8.039D+02,8.587D+02,9.534D+02,1.055D+03,1.182D+03,1.298D+03/
      DATA (SIGMA(J, 12),J=1,6) /
     &8.515D+02,8.957D+02,9.869D+02,1.122D+03,1.253D+03,1.366D+03/
      DATA (SIGMA(J, 13),J=1,6) /
     &8.769D+02,9.100D+02,1.018D+03,1.119D+03,1.252D+03,1.341D+03/
      DATA (SIGMA(J, 14),J=1,6) /
     &9.058D+02,9.532D+02,1.057D+03,1.171D+03,1.302D+03,1.391D+03/
      DATA (SIGMA(J, 15),J=1,6) /
     &9.555D+02,9.799D+02,1.098D+03,1.201D+03,1.342D+03,1.444D+03/
      DATA (SIGMA(J, 16),J=1,6) /
     &1.009D+03,1.058D+03,1.149D+03,1.290D+03,1.414D+03,1.520D+03/
      DATA (SIGMA(J, 17),J=1,6) /
     &9.907D+02,1.045D+03,1.166D+03,1.290D+03,1.384D+03,1.516D+03/
      DATA (SIGMA(J, 18),J=1,6) /
     &1.036D+03,1.121D+03,1.198D+03,1.328D+03,1.470D+03,1.592D+03/
      DATA (SIGMA(J, 19),J=1,6) /
     &1.083D+03,1.162D+03,1.250D+03,1.371D+03,1.516D+03,1.661D+03/
      DATA (SIGMA(J, 20),J=1,6) /
     &1.146D+03,1.215D+03,1.295D+03,1.443D+03,1.544D+03,1.744D+03/
      DATA (SIGMA(J, 21),J=1,6) /
     &1.158D+03,1.234D+03,1.292D+03,1.467D+03,1.618D+03,1.750D+03/
      DATA (SIGMA(J, 22),J=1,6) /
     &1.153D+03,1.205D+03,1.329D+03,1.451D+03,1.596D+03,1.734D+03/
      DATA (SIGMA(J, 23),J=1,6) /
     &1.210D+03,1.274D+03,1.356D+03,1.493D+03,1.655D+03,1.803D+03/
      DATA (SIGMA(J, 24),J=1,6) /
     &1.212D+03,1.273D+03,1.398D+03,1.489D+03,1.641D+03,1.800D+03/
      DATA (SIGMA(J, 25),J=1,6) /
     &1.236D+03,1.315D+03,1.423D+03,1.561D+03,1.669D+03,1.855D+03/
      DATA (SIGMA(J, 26),J=1,6) /
     &1.279D+03,1.345D+03,1.431D+03,1.595D+03,1.734D+03,1.889D+03/
      DATA (SIGMA(J, 27),J=1,6) /
     &1.228D+03,1.304D+03,1.438D+03,1.546D+03,1.714D+03,1.836D+03/
      DATA (SIGMA(J, 28),J=1,6) /
     &1.289D+03,1.370D+03,1.451D+03,1.597D+03,1.754D+03,1.913D+03/
      DATA (SIGMA(J, 29),J=1,6) /
     &1.411D+03,1.469D+03,1.613D+03,1.777D+03,1.910D+03,2.075D+03/
      DATA (SIGMA(J, 30),J=1,6) /
     &1.347D+03,1.401D+03,1.498D+03,1.642D+03,1.816D+03,1.975D+03/
      DATA (SIGMA(J, 31),J=1,6) /
     &1.359D+03,1.448D+03,1.551D+03,1.694D+03,1.858D+03,2.007D+03/
      DATA (SIGMA(J, 32),J=1,6) /
     &1.358D+03,1.460D+03,1.559D+03,1.698D+03,1.842D+03,1.974D+03/
      DATA (SIGMA(J, 33),J=1,6) /
     &1.418D+03,1.448D+03,1.578D+03,1.727D+03,1.872D+03,2.047D+03/
      DATA (SIGMA(J, 34),J=1,6) /
     &1.433D+03,1.466D+03,1.605D+03,1.738D+03,1.892D+03,2.019D+03/
      DATA (SIGMA(J, 35),J=1,6) /
     &1.430D+03,1.511D+03,1.602D+03,1.752D+03,1.935D+03,2.060D+03/
      DATA (SIGMA(J, 36),J=1,6) /
     &1.462D+03,1.499D+03,1.653D+03,1.805D+03,1.920D+03,2.057D+03/
      DATA (SIGMA(J, 37),J=1,6) /
     &1.470D+03,1.520D+03,1.656D+03,1.818D+03,1.946D+03,2.131D+03/
      DATA (SIGMA(J, 38),J=1,6) /
     &1.470D+03,1.542D+03,1.691D+03,1.800D+03,1.968D+03,2.133D+03/
      DATA (SIGMA(J, 39),J=1,6) /
     &1.495D+03,1.588D+03,1.676D+03,1.834D+03,1.969D+03,2.163D+03/
      DATA (SIGMA(J, 40),J=1,6) /
     &1.525D+03,1.551D+03,1.722D+03,1.833D+03,2.020D+03,2.192D+03/
      DATA (SIGMA(J, 41),J=1,6) /
     &1.526D+03,1.615D+03,1.709D+03,1.899D+03,2.040D+03,2.181D+03/
      DATA (SIGMA(J, 42),J=1,6) /
     &1.510D+03,1.567D+03,1.716D+03,1.892D+03,2.056D+03,2.197D+03/
      DATA (SIGMA(J, 43),J=1,6) /
     &1.557D+03,1.658D+03,1.776D+03,1.898D+03,2.092D+03,2.200D+03/
      DATA (SIGMA(J, 44),J=1,6) /
     &1.556D+03,1.645D+03,1.752D+03,1.920D+03,2.091D+03,2.243D+03/
      DATA (SIGMA(J, 45),J=1,6) /
     &1.583D+03,1.663D+03,1.798D+03,1.940D+03,2.051D+03,2.263D+03/
      DATA (SIGMA(J, 46),J=1,6) /
     &1.599D+03,1.642D+03,1.799D+03,1.941D+03,2.107D+03,2.268D+03/
      DATA (SIGMA(J, 47),J=1,6) /
     &1.611D+03,1.692D+03,1.811D+03,1.956D+03,2.107D+03,2.264D+03/
      DATA (SIGMA(J, 48),J=1,6) /
     &1.625D+03,1.706D+03,1.819D+03,1.986D+03,2.139D+03,2.354D+03/
      DATA (SIGMA(J, 49),J=1,6) /
     &1.666D+03,1.737D+03,1.854D+03,1.971D+03,2.160D+03,2.318D+03/
      DATA (SIGMA(J, 50),J=1,6) /
     &1.648D+03,1.747D+03,1.856D+03,2.023D+03,2.181D+03,2.352D+03/
      DATA (SIGMA(J, 51),J=1,6) /
     &1.653D+03,1.763D+03,1.868D+03,2.015D+03,2.203D+03,2.386D+03/
      DATA (SIGMA(J, 52),J=1,6) /
     &1.690D+03,1.720D+03,1.902D+03,2.027D+03,2.189D+03,2.357D+03/
      DATA (SIGMA(J, 53),J=1,6) /
     &1.690D+03,1.750D+03,1.921D+03,2.059D+03,2.208D+03,2.417D+03/
      DATA (SIGMA(J, 54),J=1,6) /
     &1.705D+03,1.781D+03,1.911D+03,2.073D+03,2.242D+03,2.411D+03/
      DATA (SIGMA(J, 55),J=1,6) /
     &1.714D+03,1.806D+03,1.896D+03,2.100D+03,2.253D+03,2.411D+03/
      DATA (SIGMA(J, 56),J=1,6) /
     &1.774D+03,1.813D+03,1.954D+03,2.098D+03,2.280D+03,2.482D+03/
 
      DATA (SIGQE(J, 2),J=1,6) /
     &4.141D+01,3.708D+01,5.428D+01,8.696D+01,1.403D+02,1.885D+02/
      DATA (SIGQE(J, 3),J=1,6) /
     &4.357D+01,3.894D+01,5.177D+01,9.675D+01,1.447D+02,2.029D+02/
      DATA (SIGQE(J, 4),J=1,6) /
     &4.123D+01,3.933D+01,6.070D+01,9.482D+01,1.474D+02,2.023D+02/
      DATA (SIGQE(J, 5),J=1,6) /
     &4.681D+01,4.287D+01,6.381D+01,1.050D+02,1.519D+02,2.198D+02/
      DATA (SIGQE(J, 6),J=1,6) /
     &5.407D+01,5.195D+01,6.723D+01,1.108D+02,1.750D+02,2.368D+02/
      DATA (SIGQE(J, 7),J=1,6) /
     &4.975D+01,4.936D+01,6.880D+01,1.162D+02,1.689D+02,2.329D+02/
      DATA (SIGQE(J, 8),J=1,6) /
     &5.361D+01,5.027D+01,6.858D+01,1.177D+02,1.759D+02,2.412D+02/
      DATA (SIGQE(J, 9),J=1,6) /
     &4.980D+01,5.063D+01,7.210D+01,1.196D+02,1.806D+02,2.299D+02/
      DATA (SIGQE(J, 10),J=1,6) /
     &5.170D+01,5.070D+01,7.105D+01,1.182D+02,1.679D+02,2.411D+02/
      DATA (SIGQE(J, 11),J=1,6) /
     &4.950D+01,4.950D+01,7.286D+01,1.137D+02,1.769D+02,2.477D+02/
      DATA (SIGQE(J, 12),J=1,6) /
     &5.262D+01,5.133D+01,7.110D+01,1.204D+02,1.789D+02,2.501D+02/
      DATA (SIGQE(J, 13),J=1,6) /
     &5.320D+01,5.378D+01,6.847D+01,1.200D+02,1.805D+02,2.442D+02/
      DATA (SIGQE(J, 14),J=1,6) /
     &5.638D+01,5.271D+01,6.985D+01,1.209D+02,1.867D+02,2.610D+02/
      DATA (SIGQE(J, 15),J=1,6) /
     &5.294D+01,5.353D+01,7.435D+01,1.211D+02,1.899D+02,2.612D+02/
      DATA (SIGQE(J, 16),J=1,6) /
     &5.668D+01,5.254D+01,7.557D+01,1.269D+02,1.917D+02,2.707D+02/
      DATA (SIGQE(J, 17),J=1,6) /
     &5.456D+01,5.721D+01,7.481D+01,1.208D+02,1.859D+02,2.658D+02/
      DATA (SIGQE(J, 18),J=1,6) /
     &5.901D+01,5.382D+01,7.591D+01,1.246D+02,1.872D+02,2.874D+02/
      DATA (SIGQE(J, 19),J=1,6) /
     &6.328D+01,6.116D+01,8.451D+01,1.318D+02,2.088D+02,2.749D+02/
      DATA (SIGQE(J, 20),J=1,6) /
     &5.779D+01,5.924D+01,8.382D+01,1.370D+02,2.062D+02,2.837D+02/
      DATA (SIGQE(J, 21),J=1,6) /
     &7.155D+01,5.732D+01,8.231D+01,1.363D+02,2.047D+02,2.820D+02/
      DATA (SIGQE(J, 22),J=1,6) /
     &6.699D+01,5.651D+01,8.511D+01,1.477D+02,2.031D+02,2.921D+02/
      DATA (SIGQE(J, 23),J=1,6) /
     &6.179D+01,6.269D+01,9.395D+01,1.437D+02,2.195D+02,2.964D+02/
      DATA (SIGQE(J, 24),J=1,6) /
     &6.784D+01,6.028D+01,8.622D+01,1.279D+02,2.214D+02,2.867D+02/
      DATA (SIGQE(J, 25),J=1,6) /
     &6.589D+01,5.795D+01,8.890D+01,1.385D+02,2.055D+02,2.988D+02/
      DATA (SIGQE(J, 26),J=1,6) /
     &6.364D+01,6.325D+01,8.942D+01,1.421D+02,2.128D+02,3.083D+02/
      DATA (SIGQE(J, 27),J=1,6) /
     &6.449D+01,6.664D+01,8.986D+01,1.453D+02,2.140D+02,2.932D+02/
      DATA (SIGQE(J, 28),J=1,6) /
     &7.284D+01,6.139D+01,8.867D+01,1.425D+02,2.179D+02,2.978D+02/
      DATA (SIGQE(J, 29),J=1,6) /
     &7.221D+01,7.085D+01,9.079D+01,1.482D+02,2.277D+02,2.913D+02/
      DATA (SIGQE(J, 30),J=1,6) /
     &6.928D+01,6.294D+01,8.935D+01,1.463D+02,2.265D+02,2.834D+02/
      DATA (SIGQE(J, 31),J=1,6) /
     &6.611D+01,6.586D+01,9.133D+01,1.461D+02,2.201D+02,2.959D+02/
      DATA (SIGQE(J, 32),J=1,6) /
     &6.401D+01,6.177D+01,8.971D+01,1.480D+02,2.155D+02,3.152D+02/
      DATA (SIGQE(J, 33),J=1,6) /
     &7.057D+01,6.918D+01,8.410D+01,1.465D+02,2.288D+02,3.088D+02/
      DATA (SIGQE(J, 34),J=1,6) /
     &6.453D+01,7.020D+01,9.272D+01,1.517D+02,2.189D+02,2.999D+02/
      DATA (SIGQE(J, 35),J=1,6) /
     &6.741D+01,6.295D+01,9.323D+01,1.536D+02,2.190D+02,2.930D+02/
      DATA (SIGQE(J, 36),J=1,6) /
     &6.807D+01,7.046D+01,1.025D+02,1.565D+02,2.315D+02,3.090D+02/
      DATA (SIGQE(J, 37),J=1,6) /
     &8.082D+01,6.565D+01,9.160D+01,1.572D+02,2.229D+02,3.125D+02/
      DATA (SIGQE(J, 38),J=1,6) /
     &6.494D+01,6.964D+01,9.089D+01,1.653D+02,2.336D+02,3.120D+02/
      DATA (SIGQE(J, 39),J=1,6) /
     &6.833D+01,6.860D+01,8.933D+01,1.601D+02,2.261D+02,3.167D+02/
      DATA (SIGQE(J, 40),J=1,6) /
     &7.021D+01,6.866D+01,8.437D+01,1.588D+02,2.249D+02,2.941D+02/
      DATA (SIGQE(J, 41),J=1,6) /
     &7.122D+01,6.205D+01,9.545D+01,1.582D+02,2.335D+02,3.395D+02/
      DATA (SIGQE(J, 42),J=1,6) /
     &7.265D+01,6.936D+01,9.486D+01,1.505D+02,2.379D+02,3.248D+02/
      DATA (SIGQE(J, 43),J=1,6) /
     &7.048D+01,7.539D+01,9.192D+01,1.566D+02,2.532D+02,3.182D+02/
      DATA (SIGQE(J, 44),J=1,6) /
     &6.650D+01,7.139D+01,9.862D+01,1.602D+02,2.289D+02,3.077D+02/
      DATA (SIGQE(J, 45),J=1,6) /
     &7.511D+01,6.893D+01,9.245D+01,1.641D+02,2.519D+02,3.381D+02/
      DATA (SIGQE(J, 46),J=1,6) /
     &6.437D+01,6.894D+01,8.697D+01,1.544D+02,2.391D+02,3.213D+02/
      DATA (SIGQE(J, 47),J=1,6) /
     &7.980D+01,6.958D+01,1.022D+02,1.609D+02,2.408D+02,3.246D+02/
      DATA (SIGQE(J, 48),J=1,6) /
     &7.265D+01,7.313D+01,8.989D+01,1.578D+02,2.387D+02,3.235D+02/
      DATA (SIGQE(J, 49),J=1,6) /
     &6.959D+01,6.337D+01,9.084D+01,1.656D+02,2.331D+02,3.226D+02/
      DATA (SIGQE(J, 50),J=1,6) /
     &7.371D+01,6.807D+01,9.726D+01,1.535D+02,2.445D+02,3.189D+02/
      DATA (SIGQE(J, 51),J=1,6) /
     &7.882D+01,6.680D+01,9.377D+01,1.629D+02,2.448D+02,3.297D+02/
      DATA (SIGQE(J, 52),J=1,6) /
     &7.223D+01,6.794D+01,9.925D+01,1.738D+02,2.446D+02,3.162D+02/
      DATA (SIGQE(J, 53),J=1,6) /
     &7.703D+01,6.971D+01,9.601D+01,1.595D+02,2.484D+02,3.265D+02/
      DATA (SIGQE(J, 54),J=1,6) /
     &7.549D+01,7.459D+01,8.984D+01,1.645D+02,2.348D+02,3.201D+02/
      DATA (SIGQE(J, 55),J=1,6) /
     &7.891D+01,6.840D+01,1.017D+02,1.698D+02,2.501D+02,3.429D+02/
      DATA (SIGQE(J, 56),J=1,6) /
     &7.545D+01,6.673D+01,1.057D+02,1.684D+02,2.424D+02,3.181D+02/

      ASQS = 0.5D0*LOG10(1.876D+03*E0)
      JE = MIN(INT((ASQS-AMIN)/DA)+1,NE-2)
      DO JA=2,IA
         ABEAM = DBLE(JA)
         S1 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
     +                   SIGMA(JE,JA),SIGMA(JE+1,JA),SIGMA(JE+2,JA))
         S2 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
     +                   SIGQE(JE,JA),SIGQE(JE+1,JA),SIGQE(JE+2,JA))
         SSIGNUC(JA) = S1 + S2
         ALNUC(JA) = ATARGET/(AVOG*SSIGNUC(JA))
      ENDDO
      ALNUC(1) = FPNI(E0, 13)
      SSIGNUC(1) = ATARGET/(AVOG*ALNUC(1))

      RETURN
      END


C*=======================================================================
C.  General utilities
C=======================================================================

      FUNCTION QUAD_INT (R,X0,X1,X2,V0,V1,V2)

C-----------------------------------------------------------------------
C...Quadratic interpolation
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      SAVE

      R0=R-X0
      R1=R-X1
      R2=R-X2
      S0=X0-X1
      S1=X0-X2
      S2=X1-X2
      QUAD_INT = V0*R1*R2/(S0*S1)-V1*R0*R2/(S0*S2)+V2*R0*R1/(S1*S2)
      RETURN
      END
C=======================================================================

      FUNCTION GAUSS (FUN, A,B)

C-----------------------------------------------------------------------
C...Returns the  8 points Gauss-Legendre integral
C.  of function FUN from A to B
C...........................................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      DIMENSION X(8), W(8)
      SAVE
      DATA X/.0950125098D0, .2816035507D0, .4580167776D0, .6178762444D0,
     1       .7554044083D0, .8656312023D0, .9445750230D0, .9894009349D0/
      DATA W/.1894506104D0, .1826034150D0, .1691565193D0, .1495959888D0,
     1       .1246289712D0, .0951585116D0, .0622535239D0, .0271524594D0/

      XM = 0.5D0*(B+A)
      XR = 0.5D0*(B-A)
      SS = 0.D0
      DO J=1,8
        DX = XR*X(J)
        SS = SS + W(J) * (FUN(XM+DX) + FUN(XM-DX))
      ENDDO
      GAUSS = XR*SS
      RETURN
      END
C=======================================================================

      SUBROUTINE INVERT_ARRAY (yy, xmin, dx, n, xnew, ymin, dy)

C-----------------------------------------------------------------------
C..    This subroutine receives one   array
C      of n y values in input yy(1:n)
C      that correspond to  equispaced values of x_j = xmin + dx*(j-1)
C
C      and "reverse" the array returning an array of  x values
C      xnew (1:n) that  corresponds to equispaced values of y
C      The relation is assumed monotonous but can be 
C      increasing or decreasing
C..............................................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      dimension  yy(n), xnew (n)
      SAVE

      ymin = yy(1)
      ymax = yy(n)
      dy = (ymax - ymin)/float(n-1)
      xnew (1) = xmin
      xnew (n) = xmin + dx*float(n-1)
      k0 = 1
      do j=2,n-1
         y = ymin + float(j-1)*dy 
         do k=k0,n
            if((yy(k) .gt. y) .eqv. (yy(n) .gt. yy(1))) goto 100
         enddo
100      y2 = yy(k)
         y1 = yy(k-1)
         k0 = k-1
         x1 = xmin + dx*float(k-2)
         x2 = x1+dx
         xnew (j)  = x1 + dx* (y-y1)/(y2-y1)
      enddo
      return
      end
C->
C=======================================================================

      SUBROUTINE SINCO(S,C)

C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

      F = TWOPI*S_RNDM(0)
      C = COS (F)
      S = SIN (F)
      RETURN
      END

C***********************************************************************
C.  Cross sections for cascade calculations (FPNI)
C=======================================================================
      
      SUBROUTINE SIGMA_PP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO) 

C-----------------------------------------------------------------------
C...p-p cross sections
C.
C.  this routine serves the purpose to calculate cascades with different 
C.  cross sections
C.
C. INPUT: E0 = Laboratory Energy  (TeV)
C. 
C. OUTPUT: SIGT = total cross section
C.         SIGEL = elastic cross section
C.         SIGINEL = inelastic cross section
C.         SLOPE = slope of elastic scattering (GeV**-2)
C.         RHO = Imaginary/Real part of forward elastic amplitude
C.   
C.  (old cross section tables end at 10^6 GeV)
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      DIMENSION SSIG0(51)
      DIMENSION SIGDIF(3)
      COMMON /CSPA/ ICSPA2(3)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

C...p-p inelastic cross sections (mbarn)
      DATA (SSIG0(J),J=1,51) /
     +      32.05D0,  32.06D0,  32.08D0,  32.13D0,  32.22D0,  32.36D0,
     +      32.56D0,  32.85D0,  33.24D0,  33.75D0,  34.37D0,  35.14D0,
     +      36.05D0,  37.12D0,  38.37D0,  39.78D0,  41.36D0,  43.13D0,
     +      45.07D0,  47.18D0,  49.47D0,  51.91D0,  54.54D0,  57.28D0,
     +      60.15D0,  63.15D0,  66.28D0,  69.48D0,  72.80D0,  76.22D0,
     +      79.71D0,  83.27D0,  86.87D0,  90.55D0,  94.26D0,  98.05D0,
     +     101.89D0, 105.75D0, 109.71D0, 113.65D0, 117.60D0, 121.55D0,
     +     125.53D0, 129.56D0, 133.60D0, 137.70D0, 141.77D0, 145.84D0,
     +     149.92D0, 154.02D0, 158.15D0/

      ICSPA = ICSPA2(1)

      SQS = SQRT(2000.D0*0.938D0*E0)

*  pre-LHC SIBYLL2.1 model
      
      IF(ICSPA.EQ.-2) THEN

         CALL SIB_SIGMA_EXT(3,SQS,SIGT,SIGEL,SIGINEL,SLOPE,RHO)      

*  old standard NUCLIB/SIBYLL model

      ELSE IF(ICSPA.EQ.-1) THEN

        AL = LOG10(SQS)
        if(AL.le.1.D0) then
          SIGINEL = SSIG0(1)
        else
          J1 = INT((AL - 1.D0)*10.D0) + 1
          J1 = min(J1,50)
          T = (AL-1.D0)*10.D0 - DBLE(J1-1)
          SIGINEL = SSIG0(J1)*(1.D0-T) + SSIG0(J1+1)*T
        endif
        CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,SIGEL1,SIGEL2)
        R = SIGEL1/SIGT1
        RHO = RHO1
        SIGT  = SIGINEL/(1.D0-R)
        SIGEL = SIGINEL*R/(1.D0-R)
        SLOPE = SIGT**2/(SIGEL * 16.D0*PI) * (1.D0+RHO1**2) /CMBARN

*  cross section as calculated in SIBYLL

      ELSE IF(ICSPA.EQ.0) THEN

        CALL SIB_SIGMA_HP(1,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)

*  Donnachie-Landshoff  (sig-tot)

      ELSE IF(ICSPA.EQ.1) THEN

        CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,
     +             SIGEL1,SIGEL2)
        R = SIGEL1/SIGT1
        RHO = RHO1

        DELDL = 0.0808D0
        EPSDL = -0.4525D0
        S = SQS*SQS
        SIGT = 21.7D0*S**DELDL+56.08D0*S**EPSDL
        SIGEL = R*SIGT
        SIGINEL = SIGT-SIGEL
        SLOPE = SIGT**2/(SIGEL * 16.D0*PI) * (1.D0+RHO**2) /CMBARN

*  Donnachie-Landshoff (sig-tot and sig-el)

      ELSE IF(ICSPA.EQ.2) THEN

        DELDL = 0.0808D0
        EPSDL = -0.4525D0
        S = SQS*SQS
        SIGT = 21.7D0*S**DELDL+56.08D0*S**EPSDL
        IMODEL = 1
        IF(IMODEL.EQ.1) THEN
          ALPHAP = 0.25D0
          SLOPE = 8.5D0+2.D0*ALPHAP*LOG(S)
        ELSE IF(IMODEL.EQ.2) THEN
          ALPHAP = 0.3D0
          SLOPE = 8.D0+2.D0*ALPHAP*LOG(S)
        ENDIF
        SIGEL = SIGT**2/(16.D0*PI*SLOPE*CMBARN)
        SIGINEL = SIGT-SIGEL
        RHO = 0.D0

*  geometrical scaling with Donnachie-Landshoff sig-tot

      ELSE IF(ICSPA.EQ.3) THEN

        R = 0.17D0

        DELDL = 0.0808D0
        EPSDL = -0.4525D0
        S = SQS*SQS
        SIGT = 21.7D0*S**DELDL+56.08D0*S**EPSDL

        SIGEL = R*SIGT
        SIGINEL = SIGT-SIGEL
        SLOPE = SIGT**2/(16.D0*PI*SIGEL)/CMBARN
        RHO = 0.D0

c ICSPA=4 reserved for CONEX_EXTENSION
c      ELSE IF(ICSPA.EQ.4) THEN

*  cross section from 2014 Review of Particle Physics
        
      ELSE IF(ICSPA.EQ.5) THEN
         
c     elastic slope not included in fit
c     taking slope parameterization from sigma_pp Donnie.-Landshoff
         ALPHAP = 0.25D0
         SLOPE = 8.5D0+4.D0*ALPHAP*LOG(SQS)
         
         CALL SIG_RPP2014(1,1,SQS,SLOPE,SIGT,SIGEL,SIGINEL,RHO)

      ENDIF

      RETURN
      END

C=======================================================================

      SUBROUTINE SIGMA_PIP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO) 

C-----------------------------------------------------------------------
C...pi-p cross sections
C.
C.  this routine serves the purpose to calculate cascades with different 
C.  cross sections
C.
C. INPUT: E0 = Laboratory Energy  (TeV)
C. 
C. OUTPUT: SIGT = total cross section
C.         SIGEL = elastic cross section
C.         SIGINEL = inelastic cross section
C.         SLOPE = slope of elastic scattering (GeV**-2)
C.         RHO = Imaginary/Real part of forward elastic amplitude
C.
C.  (old cross section tables end at 10^6 GeV)
C-----------------------------------------------------------------------
Cf2py double precision,intent(in) :: e0
Cf2py double precision,intent(out) :: sigt, sigel, siginel, slope, rho
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      DIMENSION SSIG0(51)
      DIMENSION SIGDIF(3)
      COMMON /CSPA/ ICSPA2(3)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

C...pi-p inelastic cross sections (mbarn)
      DATA (SSIG0(J),J=1,51) /
     +      20.76D0,  20.78D0,  20.81D0,  20.88D0,  20.98D0,  21.13D0,
     +      21.33D0,  21.61D0,  21.96D0,  22.39D0,  22.92D0,  23.56D0,
     +      24.31D0,  25.18D0,  26.18D0,  27.32D0,  28.60D0,  30.04D0,
     +      31.64D0,  33.40D0,  35.34D0,  37.43D0,  39.72D0,  42.16D0,
     +      44.77D0,  47.56D0,  50.53D0,  53.66D0,  56.99D0,  60.50D0,
     +      64.17D0,  68.03D0,  72.05D0,  76.27D0,  80.67D0,  85.27D0,
     +      90.08D0,  95.04D0, 100.27D0, 105.65D0, 111.21D0, 116.94D0,
     +     122.87D0, 129.03D0, 135.37D0, 141.93D0, 148.62D0, 155.49D0,
     +     162.48D0, 169.60D0, 176.94D0/

      ICSPA = ICSPA2(2)

      SQS = SQRT(2000.D0*0.938D0*E0)
      
*  pre-LHC SIBYLL2.1 model
      
      IF(ICSPA.EQ.-2) THEN

         CALL SIB_SIGMA_EXT(2,SQS,SIGT,SIGEL,SIGINEL,SLOPE,RHO)
      
*  old standard NUCLIB/SIBYLL model

      ELSE IF(ICSPA.EQ.-1) THEN

        AL = LOG10(SQS)
        if(AL.le.1.D0) then
          SIGINEL = SSIG0(1)
        else
          J1 = INT((AL - 1.D0)*10.D0) + 1
          J1 = min(J1,50)
          T = (AL-1.D0)*10.D0 - DBLE(J1-1)
          SIGINEL = SSIG0(J1)*(1.D0-T) + SSIG0(J1+1)*T
        endif
        CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,SIGEL1,SIGEL2)
        R = SIGEL1/SIGT1
        RHO = RHO1
        SIGT  = SIGINEL/(1.D0-R)
        SIGEL = SIGINEL*R/(1.D0-R)
        SLOPE = SIGT**2/(SIGEL * 16.D0*PI) * (1.D0+RHO1**2) /CMBARN

*  cross section as calculated in SIBYLL

      ELSE IF(ICSPA.EQ.0) THEN

        CALL SIB_SIGMA_HP(2,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)

*  Donnachie-Landshoff  (sig-tot)

      ELSE IF(ICSPA.EQ.1) THEN

        CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,
     +             SIGEL1,SIGEL2)
        R = SIGEL1/SIGT1
        RHO = RHO1

        DELDL = 0.0808D0
        EPSDL = -0.4525D0
        S = SQS*SQS
        SIGT = 13.63D0*S**DELDL+(36.02D0+27.56D0)/2.D0*S**EPSDL
        SIGEL = R*SIGT
        SIGINEL = SIGT-SIGEL
        SLOPE = SIGT**2/(SIGEL * 16.D0*PI) * (1.D0+RHO**2) /CMBARN

*  Donnachie-Landshoff (sig-tot and sig-el)

      ELSE IF(ICSPA.EQ.2) THEN

        DELDL = 0.0808D0
        EPSDL = -0.4525D0
        S = SQS*SQS
        SIGT = 13.63D0*S**DELDL+(36.02D0+27.56D0)/2.D0*S**EPSDL
        IMODEL = 1
        IF(IMODEL.EQ.1) THEN
          ALPHAP = 0.25D0
          SLOPE = 8.5D0+2.D0*ALPHAP*LOG(S)
        ELSE IF(IMODEL.EQ.2) THEN
          ALPHAP = 0.3D0
          SLOPE = 8.D0+2.D0*ALPHAP*LOG(S)
        ENDIF
        SIGEL = SIGT**2/(16.D0*PI*SLOPE*CMBARN)
        SIGINEL = SIGT-SIGEL
        RHO = 0.

*  geometrical scaling with Donnachie-Landshoff sig-tot

      ELSE IF(ICSPA.EQ.3) THEN

        R = 0.17D0

        DELDL = 0.0808D0
        EPSDL = -0.4525D0
        S = SQS*SQS
        SIGT = 13.63D0*S**DELDL+(36.02D0+27.56D0)/2.D0*S**EPSDL

        SIGEL = R*SIGT
        SIGINEL = SIGT-SIGEL
        SLOPE = SIGT**2/(16.D0*PI*SIGEL)/CMBARN
        RHO = 0.D0

c ICSPA=4 reserved for CONEX_EXTENSION
c      ELSE IF(ICSPA.EQ.4) THEN

*  cross section from 2014 Review of Particle Physics
        
      ELSE IF(ICSPA.EQ.5) THEN
         
c     elastic slope not included in fit
c     taking slope parameterization from sigma_pp Donnie.-Landshoff
         ALPHAP = 0.25D0
         SLOPE = 8.5D0+4.D0*ALPHAP*LOG(SQS)
         
         CALL SIG_RPP2014(2,1,SQS,SLOPE,SIGT,SIGEL,SIGINEL,RHO)
         
      ENDIF


      RETURN
      END

C=======================================================================

      SUBROUTINE SIGMA_KP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO) 

C-----------------------------------------------------------------------
C...K-p cross sections
C.
C.  this routine serves the purpose to calculate cascades with different 
C.  cross sections
C.
C.  if old cross sections are selected then sigma_pi = sigma_K
C.
C. INPUT: E0 = Laboratory Energy  (TeV)
C. 
C. OUTPUT: SIGT = total cross section
C.         SIGEL = elastic cross section
C.         SIGINEL = inelastic cross section
C.         SLOPE = slope of elastic scattering (GeV**-2)
C.         RHO = Imaginary/Real part of forward elastic amplitude
C.
C.  (old cross section tables end at 10^6 GeV)
C-----------------------------------------------------------------------
Cf2py double precision,intent(in) :: e0
Cf2py double precision,intent(out) :: sigt, sigel, siginel, slope, rho
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      DIMENSION SSIG0(51)
      DIMENSION SIGDIF(3)
      COMMON /CSPA/ ICSPA2(3)
      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN
      SAVE

C...pi-p inelastic cross sections (mbarn)
      DATA (SSIG0(J),J=1,51) /
     +      20.76D0,  20.78D0,  20.81D0,  20.88D0,  20.98D0,  21.13D0,
     +      21.33D0,  21.61D0,  21.96D0,  22.39D0,  22.92D0,  23.56D0,
     +      24.31D0,  25.18D0,  26.18D0,  27.32D0,  28.60D0,  30.04D0,
     +      31.64D0,  33.40D0,  35.34D0,  37.43D0,  39.72D0,  42.16D0,
     +      44.77D0,  47.56D0,  50.53D0,  53.66D0,  56.99D0,  60.50D0,
     +      64.17D0,  68.03D0,  72.05D0,  76.27D0,  80.67D0,  85.27D0,
     +      90.08D0,  95.04D0, 100.27D0, 105.65D0, 111.21D0, 116.94D0,
     +     122.87D0, 129.03D0, 135.37D0, 141.93D0, 148.62D0, 155.49D0,
     +     162.48D0, 169.60D0, 176.94D0/

      ICSPA = ICSPA2(3)
      
      SQS = SQRT(2000.D0*0.938D0*E0)      

*  pre-LHC SIBYLL2.1 model
      
      IF(ICSPA.EQ.-2) THEN

         CALL SIB_SIGMA_EXT(3,SQS,SIGT,SIGEL,SIGINEL,SLOPE,RHO)
      
*  old standard NUCLIB/SIBYLL model

      ELSE IF(ICSPA.EQ.-1) THEN

        AL = LOG10(SQS)
        if(AL.le.1.D0) then
          SIGINEL = SSIG0(1)
        else
          J1 = INT((AL - 1.D0)*10.D0) + 1
          J1 = min(J1,50)
          T = (AL-1.D0)*10.D0 - DBLE(J1-1)
          SIGINEL = SSIG0(J1)*(1.D0-T) + SSIG0(J1+1)*T
        endif
        CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,SIGEL1,SIGEL2)
        R = SIGEL1/SIGT1
        RHO = RHO1
        SIGT  = SIGINEL/(1.D0-R)
        SIGEL = SIGINEL*R/(1.D0-R)
        SLOPE = SIGT**2/(SIGEL * 16.D0*PI) * (1.D0+RHO1**2) /CMBARN

*  cross section as calculated in SIBYLL

      ELSE IF(ICSPA.EQ.0) THEN

        CALL SIB_SIGMA_HP(3,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)

*  Donnachie-Landshoff  (sig-tot)

      ELSE IF(ICSPA.EQ.1) THEN

        CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,
     +             SIGEL1,SIGEL2)
        R = SIGEL1/SIGT1
        RHO = RHO1

        DELDL = 0.0808D0
        EPSDL = -0.4525D0
        S = SQS*SQS
        SIGT = 11.82D0*S**DELDL+(26.36D0+ 8.15D0)/2.D0*S**EPSDL
        SIGEL = R*SIGT
        SIGINEL = SIGT-SIGEL
        SLOPE = SIGT**2/(SIGEL * 16.D0*PI) * (1.D0+RHO**2) /CMBARN

*  Donnachie-Landshoff (sig-tot and sig-el)

      ELSE IF(ICSPA.EQ.2) THEN

        DELDL = 0.0808D0
        EPSDL = -0.4525D0
        S = SQS*SQS
        SIGT = 11.82D0*S**DELDL+(26.36D0+ 8.15D0)/2.D0*S**EPSDL
        IMODEL = 1
        IF(IMODEL.EQ.1) THEN
          ALPHAP = 0.25D0
          SLOPE = 8.5D0+2.D0*ALPHAP*LOG(S)
        ELSE IF(IMODEL.EQ.2) THEN
          ALPHAP = 0.3D0
          SLOPE = 8.D0+2.D0*ALPHAP*LOG(S)
        ENDIF
        SIGEL = SIGT**2/(16.D0*PI*SLOPE*CMBARN)
        SIGINEL = SIGT-SIGEL
        RHO = 0.D0

*  geometrical scaling with Donnachie-Landshoff sig-tot

      ELSE IF(ICSPA.EQ.3) THEN

        R = 0.17D0

        DELDL = 0.0808D0
        EPSDL = -0.4525D0
        S = SQS*SQS
        SIGT = 11.82D0*S**DELDL+(26.36D0+ 8.15D0)/2.D0*S**EPSDL

        SIGEL = R*SIGT
        SIGINEL = SIGT-SIGEL
        SLOPE = SIGT**2/(16.D0*PI*SIGEL)/CMBARN
        RHO = 0.D0
        
c ICSPA=4 reserved for CONEX_EXTENSION
c      ELSE IF(ICSPA.EQ.4) THEN


*  cross section from 2014 Review of Particle Physics
        
      ELSE IF(ICSPA.EQ.5) THEN
         
c     elastic slope not included in fit
c     taking slope parameterization from sigma_pp Donnie.-Landshoff
         ALPHAP = 0.25D0
         SLOPE = 8.5D0+4.D0*ALPHAP*LOG(SQS)
         
         CALL SIG_RPP2014(3,1,SQS,SLOPE,SIGT,SIGEL,SIGINEL,RHO)

      ENDIF

      RETURN
      END

C=======================================================================

      SUBROUTINE SIGMA_INI 

C-----------------------------------------------------------------------
C.  Initialize the cross section and interaction lengths in air
C.  cross section model can be chosen, per particle, by setting ICSPA2()
C.  default is Sibyll cross section (0,0,0)      
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON /CSAIR/ ASQSMIN, ASQSMAX, DASQS,
     &     SSIG0(61,3),SSIGA(61,3),ALINT(61,3),NSQS
      
      COMMON /CSPA/ ICSPA2(3)

      INTEGER NCALL, NDEBUG, LUN
      COMMON /S_DEBUG/ NCALL, NDEBUG, LUN

C--------------------------------------------------------------------
C     SIBYLL utility common blocks containing constants       \FR'14
C--------------------------------------------------------------------
      DOUBLE PRECISION EPS3,EPS5,EPS8,EPS10
      COMMON /SIB_EPS/ EPS3,EPS5,EPS8,EPS10

      DOUBLE PRECISION PI,TWOPI,CMBARN
      COMMON /SIB_CST/ PI,TWOPI,CMBARN

      DOUBLE PRECISION FACN
      DIMENSION FACN(3:10)
      COMMON /SIB_FAC/ FACN
      SAVE
      DATA ICSPA2 /0,0,0/
      DATA AVOG /6.0221367D-04/
      DATA ATARGET /14.514D0/            ! effective masss of air

      IF(NDEBUG.gt.0)
     &     write(lun,*) ' SIGMA_INI: using cross section model no.',
     &     (ICSPA2(i),i=1,3)

      CALL BLOCK_INI

C...Loop on c.m. energy 
      NSQS = 61
      SQSMIN = 10.D0
      SQSMAX = 1.d+07
      ASQSMIN = LOG10(SQSMIN)
      ASQSMAX = LOG10(SQSMAX)
      DASQS = (ASQSMAX-ASQSMIN)/DBLE(NSQS-1)
      DO J=1,NSQS
         ASQS = ASQSMIN + DASQS*DBLE(J-1)
         SQS = 10.D0**ASQS
         E0 = SQS*SQS/(2.D0*0.938D0) * 1.D-03       ! TeV
C...p-air
         CALL SIGMA_PP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
C     using parametrization by Goulianos for diff. cross section
c     (depends on elastic cross section)
c     used to determine coupling to intermediate resonances in Glauber calc (ALAM)
c     assumed to be universal, i.e. same coupling used for proton, pion and kaons
         CALL SIB_HADCS1(1,SQS,SIGT1,SIGEL1,SIGINEL1,SLOPE1,RHO1)
         SIGEFF = 0.68D0*(1.D0+36.D0/SQS**2)
     &        *LOG(0.6D0+0.02D0/1.5D0*SQS**2)
         SIGEFF = MAX(0.D0,SIGEFF)
         ALAM = sqrt(SIGEFF/SIGEL1)
         SSIGSD = 2.D0 * SIGEFF        
         CALL SIG_H_AIR (SIGT, SLOPE, RHO, ALAM,
     &        SSIGT, SSIGEL, SSIGQE, SIGSD, SIGQSD )
         SSIGA(J,1) = SSIGT-SSIGQE ! had-air production cross section
         SSIG0(J,1) = SIGINEL   ! had-nucleon inel. cross section
         ALINT(J,1) = 1.D0/(AVOG*SSIGA(J,1)/ATARGET) ! interaction length in air
C...pi-air
         CALL SIGMA_PIP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO) 
         CALL  SIG_H_AIR (SIGT, SLOPE, RHO, ALAM,
     &        SSIGT, SSIGEL, SSIGQE, SIGSD, SIGQSD )
         SSIGA(J,2) = SSIGT-SSIGQE
         SSIG0(J,2) = SIGINEL
         ALINT(J,2) = 1.D0/(AVOG*SSIGA(J,2)/ATARGET)
C...K-air
         CALL SIGMA_KP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO) 
         CALL  SIG_H_AIR (SIGT, SLOPE, RHO, ALAM,
     &        SSIGT, SSIGEL, SSIGQE, SIGSD, SIGQSD )
         SSIGA(J,3) = SSIGT-SSIGQE
         SSIG0(J,3) = SIGINEL
         ALINT(J,3) = 1.D0/(AVOG*SSIGA(J,3)/ATARGET)
      ENDDO

      if (ndebug .gt. 0 ) THEN
        WRITE(LUN,'(1X,A)') 
     &  ' SIGMA_INI: NUCLIB interaction lengths [g/cm**2]'
        WRITE(LUN,'(1X,A)') 
     &  '     sqs,       p-air,      pi-air,     K-air'
      DO J=1,NSQS
         SQS = 10.D0**(ASQSMIN + DASQS*DBLE(J-1))
         WRITE(LUN,'(1X,1P,4E12.3)') 
     &        SQS,ALINT(J,1),ALINT(J,2),ALINT(J,3)
        ENDDO
      endif

      RETURN
      END

C=======================================================================

      FUNCTION FPNI (E,Linp)

C-----------------------------------------------------------------------
C...This function  returns the interaction length 
C.  of an hadronic particle travelling in air
C.
C.  INPUT:   E (TeV)   particle energy
C.           Linp      particle code
C.  OUTPUT:  FPNI      (g cm-2)
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
            
      COMMON /CSAIR/ ASQSMIN, ASQSMAX, DASQS,
     &     SSIG0(61,3),SSIGA(61,3),ALINT(61,3),NSQS

      DIMENSION KK(6:14)
      SAVE
      DATA KK /3*2, 4*3, 2*1/

      SQS = SQRT(2000.D0*E*0.937D0)                        ! GeV
      AL = LOG10 (SQS)
      L = abs(Linp)
      IF (AL .LE. ASQSMIN)  THEN
         FPNI = ALINT(1,KK(L))
      ELSE
         T = (AL-ASQSMIN)/DASQS
         J = INT(T)
         J = MIN(J,NSQS-2)
         T = T-DBLE(J)
         FPNI = ((1.D0-T)*ALINT(J+1,KK(L)) + T*ALINT(J+2,KK(L)))
      ENDIF
      RETURN
      END

C=======================================================================
      
      FUNCTION FSIGHAIR (E,Linp)

C-----------------------------------------------------------------------
C...This function returns the production cross section
C.  of an hadronic particle with air calculated in NUCLIB (SIGMA_INI)     
C.
C.  INPUT:   E (TeV)   particle energy
C.           Linp      particle code
C.  OUTPUT:  SIG_PROD  (mb)
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
            
      COMMON /CSAIR/ ASQSMIN, ASQSMAX, DASQS,
     &     SSIG0(61,3),SSIGA(61,3),ALINT(61,3),NSQS

      DIMENSION KK(6:14)
      SAVE
      DATA KK /3*2, 4*3, 2*1/

      SQS = SQRT(2000.D0*E*0.937D0)                        ! GeV
      AL = LOG10 (SQS)
      L = abs(Linp)
      IF (AL .LE. ASQSMIN)  THEN
         FSIGHAIR = SSIGA(1,KK(L))
      ELSE
         T = (AL-ASQSMIN)/DASQS
         J = INT(T)
         J = MIN(J,NSQS-2)
         T = T-DBLE(J)
         FSIGHAIR = ((1.D0-T)*SSIGA(J+1,KK(L)) + T*SSIGA(J+2,KK(L)))
      ENDIF     
      RETURN
      END

C=======================================================================

      SUBROUTINE INT_LEN_INI

C-----------------------------------------------------------------------
C...Initialize the interaction lengths from NUCLIB
C-----------------------------------------------------------------------
      SAVE
      
      CALL NUC_GEOM_INI                 ! nucleus profiles
      CALL SIGMA_INI                    ! initialize cross sections

      RETURN
      END
C=======================================================================
